title: Draw an Image as a Voronoi Map

text:
# Question
Title: 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 += '<br> Using <b>'+sites.length+'</b> 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;
}
```

```
<script src="http://www.raymondhill.net/voronoi/rhill-voronoi-core.js"></script>
Line width: <input id="linewidth" type="text" size="1" value="0" />
<br>
Show sites: <input id="showSites" type="checkbox" />
<br>
Show cells: <input id="showCells" type="checkbox" checked="checked" />
<br>
<br>
Dimensions: <input id="width" type="text" size="1" value="400" /> x <input id="height" type="text" size="1" value="300" />
<br>
Paste cell specifications here
<br>
<textarea id="cells" cols="60" rows="10" oninput='drawMap()'></textarea>
<br>
<br>
<div id="output"></div>
```

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

### Related Challenges

# Answer
> 112 votes
# 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:

# Answer
> 65 votes
# 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:

# Answer
> 56 votes
# 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`

# Answer
> 47 votes
# 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:**

# Answer
> 26 votes
# 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!

# Answer
> 23 votes
# 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]}]]
   ];

```

# Answer
> 13 votes
# Python + SciPy + emcee

The algorithm I've used is the following:

1. Resize images to a smallish size (~150 pixels)
2. Make an unsharp-masked image of the maximum channel values (this helps not pick up white regions too strongly).
3. Take the absolute value.
4. Choose random points with a probability proportional to this image. This chooses points either side of discontinuities.
5. 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:

# Answer
> 4 votes
## 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:

1. 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:

2. 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.
3. Color!

## Results

`N` = 100, 500, 1000, 3000

---
Tags: popularity-contest, geometry, image-processing
---


title: CodeGolf - Ignore the noise #1

text:
# Question
Title: CodeGolf - Ignore the noise #1
### Instructions

Barry is a *not so good* backend developer in charge of a small API that should give you information about purchases made in a client's shop. However, he hasn't done a great job and your boss tells you to fix it client-side instead. Ideally you should receive comma-separated values such as `927,2,45,90` which correspond to something like `item_id`,`item_amount`,`unit_price`,`total`

In this first puzzle we only care about `item_id` and `item_amount` but we need the other fields in place to present the problem. Barry sometimes gets things mixed up and returns noise as part of the output, he also gets the order wrong, returning `noise`, `noise`, `unit_price`, `total`, `item_id`, `item_amount`.

### Your Task

You need to store in a map (or similar structure) the pairs of `item_id` and `item_amount` regardless of whether Barry returns them correctly or not and print each element in a new line in the least number of characters possible. (Storing is not mandatory, we just care about the output)

### Sample input (mixed correct and messy formats)

```
103,2,50,100
106,1,900,900
459,40,150,300,67,2
4,20,30,6000

```

In other words input will be either `a,b,x,x` or `x,x,x,x,a,b` where what we care about is `a` & `b`. You need to provide code assuming we have a variable called `G` (or any other name) with all the lines of csv.

### Sample output

```
103,2
106,1
67,2 //This one was messy
4,20

```

### Tolerance

There is a certain margin for tolerance when it comes to the answer. Answers in similar formats but adhering to the correct values, will also be accepted. Formats like `[a,b]`,`(a,b)` or `{a,b}` are valid, but `a,b` is preferred.

> Although the first puzzle had mixed acceptance due to the fact that it was easy and fast to solve, I also felt that some people liked that too. So I'm going to continue making quickgolfs for now

# Answer
# Pyth, 10 bytes

```
FNG<.<N4 2

```

Try it online: Pyth Compiler/Executor

(The first 4 chars `=G.Q` reads all the data and stores it in `G`)

### Explanation

```
FNG              for row N in G:
    .<N4            cyclic leftshift by 4
   <     2          print the first two elements

```
> 8 votes

# Answer
# Python 3, ~~123~~ ~~76~~ ~~72~~ 62 bytes

I'm a newcomer, gotta start somewhere...

```
 for l in G:
  d=l.rstrip().split(',');x=[0,4][len(d)>5];print(d[x:x+2])

```

With xnor's and DLosc's suggestions (62 bytes):

```
for l in G:d=l.rstrip().split(',');print(d[4*(len(d)>5):][:2])

```
> 7 votes

# Answer
# Perl, 24

```
perl -pe's/(\d+,){4}|,\d+,\d+$//'

```

Example:

```
perl -pe's/(\d+,){4}|,\d+,\d+$//' <in
103,2
106,1
67,2
4,20

```
> 7 votes

# Answer
# Sed, 32 bytes

```
s/,/=/4
s/.*=//
s/,/:/2
s/:.*//

```

This deletes everything up to the 4th comma if there is one. Then it deletes everything from the 2nd comma (which may have previously been 6th).

And I even managed to right-justify!
> 4 votes

# Answer
# Javascript, 44 bytes

Assuming input in variable `G`.

*Edit: I just realized, this looks a lot like nutki's perl solution.*

```
alert(G.replace(/(\d+,){4}|,\d+,\d+$/gm,''))

```

**Snippet Demo:**

```
function run(){
    G = input.value;
    /* start solution */
    alert(G.replace(/(\d+,){4}|,\d+,\d+$/gm,''))
    /* end solution */
}
```

```
<textarea id="input" cols="25" rows="6">103,2,50,100
106,1,900,900
459,40,150,300,67,2
4,20,30,6000</textarea><br />
<button id="run" onclick="run();">Run</button>
```
> 4 votes

# Answer
# Scala, 68 bytes

```
G.map(_.split(",")).map(x=>println((x++x).slice(4,6).mkString(",")))

```

-- MrBones

# Scala, ~~46~~ 70 bytes

*Edit: Print out the result*

```
for{l<-G;r=l.split(",");x=(r++r).slice(4,6).mkString(",")}{println(x)}

```

Demo
> 4 votes

# Answer
# CJam, ~~18~~ ~~15~~ 14 bytes

```
r{',/4m<2<pr}h

```

Try it online in the CJam interpreter.

### How it works

```
r      e# Read a whitespace separated token from STDIN.
{      e# Do-while loop:
  ',/  e#    Split at commas.
  4m<  e#    Rotate the resulting array 4 units to the left.
  2<   e#    Discard everything but the first 2 elements.
  p    e#    Print the array.
  r    e#    Read a whitespace separated token from STDIN.
}h     e# If the token is a nonempty string, repeat the loop.

```
> 3 votes

# Answer
# Javascript, ~~117~~ ~~112~~ ~~99~~ ~~95~~ ~~84~~ ~~76~~ 71 bytes

```
for(r of G.split('\n'))s=r.split(','),console.log(s[t=s[4]?4:0],s[t+1])

```

JS Fiddle
> 3 votes

# Answer
# POSIX shell, 30 bytes

Assuming POSIX `sed` and `cut` are present:

```
sed 's/.*/&,&/'|cut -d, -f5,6

```

Duplicates the line, leaving the data of interest ready to be extracted from fields 5 and 6.
> 3 votes

# Answer
# Pip, 12 + 1 = 13 bytes

Uses the `-l` flag, which gives the output as space-separated lists of numbers on each line.\*

```
_@[4 5]Mg^',

```

Explanation:

* The list of numbers is expected in `g`. Since `g` is normally initialized from the command-line arguments, you can also put the input there for the same results.
* `^',` splits (each string in) the list on commas.
* `M` maps a function to each element in the result.
* `_@[4 5]` defines a lambda function. The argument to the function is represented by `_`; `@[4 5]` returns a list of its 4th and 5th elements. If the list has six elements, these will be the last two. If it has four, the indices wrap around, and so the result is the first two. (E.g. `"abcd"@4 == "abcd"@0 == "a"`.)
* We now have a list of lists, like `[[103 2] [106 1] [67 2]]`, which is auto-printed at the end of the program. The `-l` flag joins each sub-list on space and then the whole thing on newline, so the result is as follows:

```
C:\>pip.py -le "_@[4 5]Mg^'," 103,2,50,100 106,1,900,900 459,40,150,300,67,2
103 2
106 1
67 2

```

<sub>* Currently. I may change this behavior, but I'm still planning to have some flag that does the same thing.</sub>
> 2 votes

# Answer
# Perl, 37

```
@c=/(\d+,\d+)/g;$_=@c>2?$c[2]:$c[0]

```

`35` characters `+2` for `-p` and `-l`. I'm not sure if there is an approach in Perl that would beat @nutki's solution, but I figured I'd post this idea I had anyway.

Run with:

```
cat input.txt | perl -ple'@c=/(\d+,\d+)/g;$_=@c>2?$c[2]:$c[0]'

```
> 1 votes

# Answer
# bash, 54

```
while IFS=, read {a..f};do echo ${e:-$a},${f:-$b};done

```
> 1 votes

# Answer
# C, 95 bytes

```
f(a,b,c,d){scanf("%d,%d,%*d,%*d,%d,%d",&a,&b,&c,&d)>2?printf("%d,%d",c,d):printf("%d,%d",a,b);}

```

The code is simple. Here is the ungolfed version:

```
f(a,b,c,d){                                    //Four int variables
    scanf("%d,%d,%*d,%*d,%d,%d",&a,&b,&c,&d)>2 //Scan input
        ?                                      //If scanf returned a value > 2
        printf("%d,%d",c,d)                    //Print last two values
        :                                      //else
        printf("%d,%d",a,b);                   //Print first two values
}

```

# Test it here
> 0 votes

# Answer
# Ruby, 53

```
G.each_line{|x|g=x.split',';p[g[4]||g[0],g[5]||g[1]]}

```

or is good.
> 0 votes

---
Tags: code-golf, string, sorting
---


title: Reverse sentences

text:
# Question
Title: Reverse sentences
The goal of this challenge is to receive an input and output that input but with sentence order reversed. Example Input:

```
Hello friend. What are you doing? I bet it is something I want to do too!

```

Example Output:

```
I bet it is something I want to do too! What are you doing? Hello friend.

```

As you can tell from the examples, your program has to deal with question marks, exclamation points, and periods. You can assume each sentence has a punctuation and than a space before the next sentence. Trailing spaces/newlines are ok, as long as it is readable. Shortest code wins.

Good luck!

EDIT: You can assume the sentences has no quotes or parentheses, but if you make your code be able to deal with both of those then you get **-5 bytes** Example output for parens/quotes:

```
"Hello, " she said. (I didn't know what she was talking about.) --> (I didn't know what she was talking about.) "Hello, " she said.

```

# Answer
# Julia, ~~45~~ 42 bytes - 5 bonus = 37

```
s->join(reverse(split(s,r"[.?!]\K "))," ")

```

This creates an anonymous function that accepts a string as input and returns the string with its sentences reversed. This handles any special character appropriately, though double quotes and dollar signs must be escaped, otherwise they aren't valid strings in Julia.

Ungolfed + explanation:

```
function f(s)
    # Get individual sentences by splitting on the spaces that
    # separate them. Spaces are identified by matching punctuation
    # then moving the position beyond that match and matching a
    # space. This is accomplished using \K.

    sentences = split(s, r"[.?!]\K ")

    # Reverse the order of the array of sentences.

    reversed_order = reverse(sentences)

    # Join the array elements into a string, separated by a space.

    join(reversed_order, " ")
end

```

Examples:

```
julia> f("Hello friend. What are you doing? I bet it is something I want to do too!")
"I bet it is something I want to do too! What are you doing? Hello friend."

julia> f("\"Hello, \" she said. (I didn't know what she was talking about.)")
"(I didn't know what she was talking about.) \"Hello, \" she said."

```

And if you don't like looking at the escaped quotes in the output:

```
julia> println(f("\"Hello, \" she said. (I didn't know what she was talking about.)"))
(I didn't know what she was talking about.) "Hello, " she said.

```

Saved 3 bytes on the regular expression thanks to Martin Büttner! Previously this used a lookbehind: `(?<=[.?!])`.
> 9 votes

# Answer
# CJam, ~~23~~ 22 bytes

I am not sure if this qualifies for the bonus or not, but here goes the solution:

```
Sq{1$".?!"-{])[}|}%]W%

```

**Code Expansion (bit outdated)**:

```
Sq+                      e# Read the input and prepend with a space
   {            }%       e# For each input character
    _".?!"&              e# Copy and check if its one of ., ? and !
           {][}&         e# If it is one of the above, wrap everything till now in an array
                         e# and start a new array to be wrapped next time we get one of those
                         e# three characters. We now have an array of strings, each having
                         e# a single sentence
                  W%     e# Reverse the ordering of these sentences
                    s(   e# Convert to string and remove the first space

```

Try it online here
> 7 votes

# Answer
# J, 35 32

It almost handles bonus input, except I have to escape single apostrophes, so I guess it doesn't count. (Also, my first submission here)

```
f=.;@|.@(]<;.2~'.?!'e.~])@,~&' '

```

Usage:

```
f 'Hello friend. What are you doing? I bet it is something I want to do too!'

```
> 5 votes

# Answer
# Perl, 27/25

```
#!perl -n
print reverse/ |[^.!?]*./g

```

Or from the command line:

```
$perl -nE'say reverse/ |[^.!?]*./g'

```
> 4 votes

# Answer
## PHP, 60

```
echo join(' ',array_reverse(preg_split('/(?<=[?!.])/',$s)));

```
> 2 votes

# Answer
# Bash + coreutils, 40 bytes

```
sed 's/\([?.!]\) */\1\n/g'|tac|tr \\n \

```

This reads from STDIN so input may me redirected in from a file, or simply piped in, e.g.:

```
$ printf 'Hello friend. What are you doing? I bet it is something I want to do too!' | sed 's/\([?.!]\) */\1\n/g'|tac|tr \\n \
I bet it is something I want to do too! What are you doing? Hello friend.
$

```
> 2 votes

# Answer
# Pip, 25 bytes

```
a.:sRV(a^@2+$ALa@*^".?!")

```

After appending a space to the input string, we find all indices of `.`, `?`, and `!`, add 2, and use the `^@` split-at operator to break the string into sentences (each with a trailing space). Reverse the list, and it's auto-printed at the end of the program. Voilà!

Example showing the stages of the main computation with input `A! B? C. D!`:

```
              ^".?!"     ["." "?" "!"]
           a@*           [[7] [4] [1 10]]
        $AL              [7 4 1 10]
      2+                 [9 6 3 12]
   a^@                   ["A! " "B? " "C. " "D! "]
RV(                 )    ["D! " "C. " "B? " "A! "]

                         D! C. B? A!

```
> 2 votes

# Answer
# Retina, ~~61~~ ~~34~~ ~~33~~ 30 bytes

*Credits to nutki for cutting this down by 24 bytes.*

```
^
#
+`(#.*[.!?]) (.+)
$2 $1
#
<empty>

```

Where `<empty>` stands for an empty line. This assumes that `#` is not part of the input, but if that's not legitimate, I could swap it out for any other character, including `"` (which I'd only need to handle for the bonus) or something unprintable. You can run the code like that in a single file if you use the `-s` flag, or you can put each line in a separate file and pass them all to Retina.

Reversing this with a single regex replacement is possible, but really cumbersome. Even with .NET balancing groups I needed something around 90 bytes, so instead I tried doing it in multiple steps.

In Retina, every pair of lines is one replacement stage, where the first line is the pattern and the second line is the replacement.

```
^
#

```

This stage simply prepares the string for further processing. It prepends a `#` as a marker. This marker indicates that everything in front of it has already been put in the right place, and everything after it still needs to be processed.

```
+`(#.*[.!?]) (.+)
$2 $1

```

This stage swaps the sentences, by repeatedly moving the last sentence in front of the `#` (which moves forwards through the string in the process). The `+`` instructs Retina to repeat this stage until the output stops changing. As an example, here is how the input `foo. bar! blah?` would be processed:

```
#foo. bar! blah?
blah? #foo. bar!
blah? bar! #foo.

```

And finally we simply remove the marker:

```
#
<empty>

```
> 2 votes

# Answer
# JavaScript (ES6) 47 ~~45~~

As it's stated, it's a simple regex exercise. In javascript:

```
// ES6 - FireFox only
F=t=>t.match(/\S[^.!?]+./g).reverse().join(' ')

// ES5 - so much longer
function G(t){return t.match(/\S[^.!?]+./g).reverse().join(' ')}

// TEST

alert(G("Hello friend. What are you doing? I bet it is something I want to do too!"))

```
> 1 votes

# Answer
# Python 2, 62

Not going to improve for the bonus, as it's probably not worth the byte cost.

```
import re
print' '.join(re.split('(?<=[?!.]).',input())[::-1])

```
> 1 votes

# Answer
# Java, 113

```
s->{String t[]=s.split("(?<=[\\.?!]) "),u="";for(int i=t.length;i-->0;)u+=t[i]+" ";return u.replaceAll(".$","");}

```
> 1 votes

# Answer
# Ruby, 48(42 without the puts) bytes

reverse\_sentence.rb

```
puts $*[0].scan(/\S[^.!?]+./).reverse.join(" ")

```

Usage:

```
ruby reverse_sentence.rb 'Hello friend. What are you doing? I bet it is something I want to do too!'

```

Output:

```
I bet it is something I want to do too! What are you doing? Hello friend.

```

Criticism more than welcome.
> 1 votes

# Answer
# Matlab (93 bytes)

```
y=[32 input('','s')];y=sortrows([cumsum(ismember(y,'?!.'),'reverse');y]',1)';disp(y(4:2:end))

```
* This assumes the input doesn't contain leading or trailing spaces
* Uses standard input and output
* Tested in Matlab 2014b
> 1 votes

# Answer
# Ruby 41

The other Ruby answers don't have enough WTF.

```
#!ruby -apF(?<=[.!?])\s
$_=$F.reverse*" "

```

This at least works in Ruby 2. If the `a` and `F` switch works in 1.8.7, I guess you could drop `$_=` to save three characters.

Reverses every line on stdin and print to stdout:

```
$ ruby foo.rb <<< "Hello. Hi. How are you? Good, you? fine, thanks."
fine, thanks. Good, you? How are you? Hi. Hello.

```
> 1 votes

# Answer
# k, 31

```
{1_,/|(0,1+&x in"?!.")_x:" ",x}

```

.

```
k){1_,/|(0,1+&x in"?!.")_x:" ",x} "Hello friend. What are you doing? I bet it is something I want to do too!"
"I bet it is something I want to do too! What are you doing? Hello friend."

```
> 1 votes

# Answer
**C# - LINQPAD - 93 - 5 = 88 bytes**

```
void Main(){string.Join(" ",Regex.Split(Console.ReadLine(),"(?<=[.?!]) ").Reverse()).Dump();}

```

**C# Console App 189 - 5 = 184 bytes**

```
using System;using System.Linq;using System.Text.RegularExpressions;class P{static void Main(){Console.WriteLine(string.Join(" ",Regex.Split(Console.ReadLine(), "(?<=[.?!]) ").Reverse()));}}

```

*regex shamelessly flogged from Alex A. :)*
> 0 votes

# Answer
# Clojure - 44 71 chars

```
(defn rs[s](reverse(re-seq #"\S.+?[.!?]"s)))

```

Improved and simplified RE, eliminated unnecessary whitespace.

Output is a sequence of the sentences in the original string, with the order of the sentences reversed:

Input: "Hello friend. What are you doing? I bet it is something I want to do too!" Output: ("I bet it is something I want to do too!" "What are you doing?" "Hello friend.")
> 0 votes

# Answer
# Ruby, 47

```
$><<gets.strip.split(/(?<=[.?!]) /).reverse*' '

```

credits to Martin Büttner, for saving some characters.
> 0 votes

# Answer
# CJam, 21 bytes

```
1q{1?_"!.?"-}%1a/W%S*

```

This works by turning spaces after `!`s, `.`s and `?`s into the number 1 (not the *character* 1 nor the character with code point 1, so the input can still contain those), splitting at 1's, reversing the order of the resulting chunks and joining by spaces.

Try it online in the CJam interpreter.

### How it works

```
1                     e# B := 1
 q                    e# Q := input()
  {         }%        e# for each C in Q (map):
   1?                 e#   C := B ? C : 1
     _"!.?"-          e#   B := string(C).strip("!.?")
              1a/     e# Q := Q.split([1])
                 W%   e# Q := reverse(Q)
                   S* e# Q := Q.join(" ")
                      e# print(Q)

```
> 0 votes

---
Tags: code-golf
---


title: Selectively murder positive integers

text:
# Question
Title: Selectively murder positive integers
**Introduction**

Arithmetic Gaol is a special facility that incarcerates positive integers. However, recently, the positive integers have been trying to escape. Therefore the wardens have decided to, um, *eliminate* some of the positive integers to send a message to the other integers. They have hired a software engineer to write a program to figure out which integers to eliminate for maximum effect.

**Input Description**

Input is given via STDIN, command line arguments, or a user input function (such as `raw_input`). You can't have it as a function argument or a preinitialized variable (e.g. this program expects input in a variable `x`).

The first line of input contains a single positive integer `n` where `8 >= n >= 3`. Following that are `n` lines containing `n` characters from the set `[1,2,3,4,5,6,7,8,9]`. Here is an example input:

```
5
22332
46351
65455
24463
65652

```

**Output Description**

The wardens would like to eliminate numbers so that the following conditions are met:

* In each row and column of the resulting grid, no number will appear twice;
* No two eliminated numbers may be adjacent horizontally or vertically;
* The surviving numbers must form an orthogonally contiguous group -- it will be possible to travel from any surviving number to any other surviving number moving only horizontally and vertically and never crossing any eliminated number.

Output the input (minus the first line), with the eliminated numbers replaced with `#`.

There may be more than one solution. If that is the case, you may output any solution.

There may also be no solution. If that is the case, output the string `no answer`.

Here is a possible output for the example input:

```
#2#3#
46351
6#4#5
24#63
#56#2

```

**Example Inputs and Outputs**

There are multiple outputs for each input, so these outputs are just examples.

Input:

```
5
46551
51565
32654
14423
43244

```

Output:

```
46#51
#156#
326#4
1#423
#324#

```

Input:

```
7
7183625
1681563
5238564
8786268
1545382
3814756
5325345

```

Output:

```
71#362#
#6815#3
5238#64
#7#62#8
154#382
3814756
#325#4#

```

Input:

```
8
21534768
75196287
68392184
96244853
44865912
76516647
89751326
43698979

```

Output:

```
21#34768
#5196287
683#21#4
9#24#853
#4865912
7#51#64#
89751326
436#8#7#

```

Input:

```
4
2222
2331
3112
1322

```

Output:

```
no answer

```

# Answer
# Ruby, ~~346~~ ~~344~~ ~~329~~ 316 bytes, sl∞∞∞∞∞∞w

This is code-golf, not code-fast, so...

```
N=/!/=~G=$*[1..-1]*?!
M=N*N
g=->i{(!G[i]||G[i]<?*||i<0||A[i])&&break
A[i]=0
[1,N+1,-1,-1-N].map{|a|g[i+a]}}
f=->w,a{A=[];g[/\d/=~G=a];/#.{#{N}}?#/!~G&&/(\d)([^!]*|.{#{N}}.{#{O=N+1}}*)\1/!~G&&A.count(0)+w==M&&N.times.map{|m|puts G[m*(1+N),N]}&&exit
(w...M).map{|j|b=a+'';b[j]=?#
w<M&&f[w+1,b]}}
f[0,G]
$><<"no answer"

```

Use it like:

```
mad_gaksha@madlab ~/Applications/Tools $ ruby -W0 c50442.rb 3 323 312 231
#23
312
231

```

The flag `-W0` is not necessary, but I suggest you either add it to disable warnings or redirect `stderr`...

Tell me if you had enough patience to run it on the examples for n=6,7,8

# Changelog

* `each``map`, thanks to @NotThatCharles
* check for adjacent deletions and same digits by two `regexp`s, similar to what @NotThatCharles suggested
* optimized reading input a little
* smaller & slower
> 0 votes

# Answer
# Ruby - ~~541...,~~ 394

The basic algorithm is a recursive depth-first search of duplicates to affirmatively select, looking through row 1, then column 1, then row 2, etc, and checking that two neighbors are not killed and that the grid is connected (that's the `break if` clause in there, and the bit that comes before it).

```
K=(0...(N=gets.to_i)*N).to_a
J=gets(p).split*''
H=->m{K&[m+1,m-1,m+N,m-N]}
Q=->k{s=[k[j=0]]
(j=s.size
s.map{|x|(s+=H[x]&k).uniq!})while s[j]
break if(K-k).any?{|m|(H[m]-k)[0]}||k!=k&s
$><<K.map{|m|[k.index(m)?J[m]:?#,m%N>N-2?"
":p]}*''|exit if !g=((0...N*2).map{|x|(k.select{|m|m.divmod(N)[x/N]==x%N}.group_by{|m|J[m]}.find{|l,c|c[1]}||[])[1]}-[p]).min
g.map{|d|Q[k-g<<d]}}
Q[K]
puts"no answer"

```

### Some neat tricks:

`if w[1]` is much shorter than `if !w.one?` and if you know there's at least one member, it's the same result.

Similarly, `[0]` is shorter than `any?` if there's no block, and `s[j]` is a cute shortcut for `j<s.size` (technically, it's more like `j.abs<s.size`)

And `y%N+(y/N).i` is much shorter than `Complex(y%N,y/N)`

Also, when there's two complicated conditionals to generate strings, it might be shorter to do `[cond1?str1a:str1b,cond2?str2a:str2b]*''` than to add all the parens or `#{}`s in.

## Ungolfing and explanation:

(This is from the 531 byte version. I've made changes. Most notably, I've since eliminated the call to product - just solve one digit per row/column at a time, and J is now just an array, indexed by integers. All coordinates are just integers.)

`H` calculates neighbors

```
def H m
  # m, like all indices, is a complex number
  #    where the real part is x and the imaginary is y
  # so neighbors are just +/-i and +/-1

  i='i'.to_c
  neighborhood = [m+1, m-1, m+i, m-i]

  # and let's just make sure to eliminate out-of-bounds cells
  K & neighborhood
end

```

`N` is the size of the grid

```
N = gets.to_i

```

`K` are the keys to the map (complex numbers)

```
# pretty self-explanatory
# a range of, e.g., if N=3, (0..8)
# mapped to (0+0i),(1+0i),(2+0i),(0+1i),(1+1i),(2+1i),...
K = (0..N**2-1).map{|y| (y%N) +(y/N).i }

```

`J` is the input map (jail)

```
# so J is [[0+0,"2"],[0+1i,"3"],....].to_h
J=K.zip($<.flat_map {|s|
  # take each input line, and...
  # remove the "\n" and then turn it into an array of chars
  s.chomp.chars
}).to_h

```

`k` are the non-killed keys

```
# starts as K

```

`Q` is the main recursive method

```
def Q k
  j=0 # j is the size of mass
  # the connected mass starts arbitrarily wherever k starts
  mass=[k[0]]
  while j < s.size # while s hasn't grown
    j = mass.size
    mass.each{|cell|
      # add all neighbors that are in k
      (mass+=H[cell] & k).uniq!
    }
  end
  # if mass != k, it's not all orthogonally connected
  is_all_connected = k!=k&mass

  # (K-k) are the killed cells
  two_neighbors_killed = (K-k).any?{|m|
    # if any neighbors of killed cells aren't in k,
    # it means it was killed, too
    (H[m]-k)[0]
  }
  # fail fast
  return if two_neighbors_killed || is_all_connected

  def u x
    x.group_by{|m|J[m]}.select{|l,c|c[1]}
  end

  rows_with_dupes = Array.new(N){|r|u[k.select{|m|m.imag==r}]}

  cols_with_dupes = Array.new(N){|r|u[k.select{|m|m.real==r}]}

  # dupes is an array of hashes
  # each hash represents one row or column.  E.g.,
  #   {
  #     "3"=>[(0+0i),(1+0i),(3+0i)],
  #     "2"=>[(2+0i),(4+0i)]
  #   }
  # means that the 0th, 1st and 3rd cells in row 0
  # all are "3", and 2nd and 4th are "2".
  # Any digits without a duplicate are rejected.
  # Any row/col without any dupes is removed here.
  dupes = (rows_with_dupes+cols_with_dupes-[{}])

  # we solve one row at a time
  first_row = dupes[0]

  if !first_row
    # no dupes => success!
    J.map{|m,v|k.member?(m)?v:?#}.each_slice(N){|s|puts s*''}
    exit
  else
    # the digit doesn't really matter
    t=first_row.values

    # cross-multiply all arrays in the row to get a
    # small search space. We choose one cell from each
    # digit grouping and drop the rest.
    t.inject(:product).map{ |*e|
      # Technically, we drop all cells, and add back the
      # chosen cells, but it's all the same.
      new_k = k-t.flatten+e.flatten

      # and then search that space, recursively
      Q[new_k]
    }
  end
end

```

The code is executed with:

```
# run with whole board
Q[K]

# if we get here, we didn't hit an exit, so we fail
puts"no answer"

```

## Changelog

**394** added @blutorange's suggestion below, and chopped out a lot more manipulation

**408** revised output once more. Also use `.min` instead of `.inject(:+)` since I'm just taking one row anyway.

**417** shorter output calculation

**421** dropped complex numbers. Just use integers. Save a bundle

**450** more input improvements

**456** input improvements

**462** incremental improvements - esp. `find`, not `select`

**475** dropped `u` and squashed the row/col dupe builder

**503** only solve one duplicate digit per row/column at a time.

**530** use `map &:pop` instead of `values`

**531** pull out the lambda that makes the dupes array

**552** oops! missed a requirement

**536** marginally improved population of dupes array (what was formerly `d`)

**541** initial
> 5 votes

# Answer
# HTML + JavaScript (*ES6*) 459

Using an HTML textarea to get multiline input.

To test, run the snippet in Firefox. Enlarge the textarea if you like, paste the input inside the textarea (beware: exact input, no extra spaces on any line), and tab away. The result is appended.

Method: a recursive Depth First Serarch. It finds non-optimal solutions for all the test cases in seconds (it's gode golf, but a valid answer *should* terminate for common test cases)

```
<textarea onchange="s=this.value,
  z=s[0],o=-~z,k=[],X='no answer',f='#',
  (R=(s,t,h={},r=[],d=0)=>(
    s.map((v,i)=>v>0&&[i%o,~(i/o)].map(p=>h[p+=v]=[...h[p]||[],i])),
    [h[i][1]&&h[i].map(p=>r[d=p]=p)for(i in h)],
    d?r.some(p=>(
      (s[p+o]!=f&s[p-o]!=f&s[p-1]!=f&s[p+1]!=f)&&
      (
        g=[...s],g[p]=f,e=[...g],n=0,
        (F=p=>e[p]>0&&[1,-1,o,-o].map(d=>F(p+d),e[p]=--n))(p<o?p+o:p-o),
        t+n==0&&!k[g]&&(k[g]=1,R(g,t-1))
      )
    )):X=s.join('')
  ))([...s.slice(2)],z*z-1),this.value+=`

`+X"></textarea>
```

**First try**

Method: a Depth First Serarch, non-recursive, using a user stack.
The function could be easily transformed in a Breadth First Search, chaging `l.pop` to `l.shift`, so using a queue instead of a stack, but it's too slow and I'm not sure it can find an optimal solution anyway.

```
S=s=>{
  z=s[0],o=-~z,s=[...s.slice(2)];
  k=[];
  l=[[z*z,s]];

  for(d=1;d&&(t=l.pop());)
  {
    [t,s]=t,--t,

    s.map((v,i)=>v>0&&[i%o,~(i/o)].map(p=>h[p+=v]=[...h[p]||[],i]),h={},r=[])

    d=0
    for(i in h)h[i][1]&&h[i].map(p=>r[d=p]=p)

    r.map(p=>
      (s[p+o]!='#'&s[p-o]!='#'&s[p-1]!='#'&s[p+1]!='#')&&
      (
        g=[...s],g[p]='#',
        f=[...g],n=0,
        (F=p=>f[p]>0&&[1,-1,o,-o].map(d=>F(p+d),f[p]=--n))(p<o?p+o:p-o),
        t+n||k[g]||(k[g]=l.push([t,g]))
      )
    )
  }
  return d?'no answer':s.join('')
}

// TEST
cases=[
 ['5\n22332\n46351\n65455\n24463\n65652','#2#3#\n46351\n6#4#5\n24#63\n#56#2']
,['5\n46551\n51565\n32654\n14423\n43244','46#51\n#156#\n326#4\n1#423\n#324#']
,['7\n7183625\n1681563\n5238564\n8786268\n1545382\n3814756\n5325345'
 ,'71#362#\n#6815#3\n5238#64\n#7#62#8\n154#382\n3814756\n#325#4#']
,['8\n21534768\n75196287\n68392184\n96244853\n44865912\n76516647\n89751326\n43698979'
 ,'21#34768\n#5196287\n683#21#4\n9#24#853\n#4865912\n7#51#64#\n89751326\n436#8#7#']
,['4\n2222\n2331\n3112\n1322','no answer']
]

  var out = ''
  var nt=0
  var exec =_=>{
    var t=cases[nt++];
    if (t) {
      var r = S(t[0]).replace(/#/g,'<b>#</b>');
      var k = t[1].replace(/#/g,'<b>#</b>');
      out += '<tr><th>Input<th>Result<th>Expected</tr>'+
      '<tr><td>'+t[0]+'</td><td>'+r+'</td><td>'+k+'</td></tr>';
      T.innerHTML = out + 'test ' + nt + ' done';
      setTimeout(exec, 100);
    }
    else
    {
      T.innerHTML = out;
    }
  }
  T.innerHTML = "...wait ..."
  setTimeout(exec, 100);
```

```
table {
  font-size: 12px;
}
td {
  font-family: monospace;
  white-space: pre;
  vertical-align: bottom;
  margin: 2x 4px;
}

b {
  color: #ccc;
}
```

```
<table id=T>
</table>
```
> 3 votes

---
Tags: code-golf, grid
---


# Question
Title: Most concise way to fold tree structures
## Background (F#)

Let there be trees:

```
type Tree<'T> = Node of 'T * Tree<'T> list

```

Now lets fold them nicely with a function called...

```
foldTree f:('a -> 'b -> 'c) -> g:('c -> 'b -> 'b) -> a:'b -> t:Tree<'a> -> 'c

```

...taking two functions `f` and `g`, an initial state `a` and of course a tree structure `t`. Similar to the well known function fold which operates on lists, this function should "merge" siblings with `g` and parents with their children with `f` resulting in an accumulated simple value.

## Example (F#)

The tree...

```
//      1
//     / \
//    2   3
//   / \   \
//  4   5   6

let t = Node (1, [Node (2, [Node (4, []); Node (5, [])]); Node (3, [Node (6, [])])])

```

...passed to `foldTree` with the operators for addition and multiplication along with the initial state `1`...

```
let result = foldTree (+) (*) 1 t
        // = (1 + ((2 + ((4 + a) * ((5 + a) * a))) * ((3 + ((6 + a) * a)) * a)))
        // = (1 + ((2 + ((4 + 1) * ((5 + 1) * 1))) * ((3 + ((6 + 1) * 1)) * 1)))

```

...should return the value `321` to `result`.

## The Challenge

In any programming language, define the function `foldTree` in the most concise way you can come up with. Fewest number of characters wins.

# Answer
## Haskell, ~~37~~ 35:

```
data T a=a:*[T a]

(f%g)z(x:*t)=x`f`foldr(g.(f%g)z)z t

```

not counting the type definition. With it, ~~54~~ **52** (shortened by using infix operator, similarly to the answer by *proud haskeller*, but in a different way).

Ungolfed:

```
data Tree a = Node a [Tree a]

foldTree :: (a -> b -> c) -> (c -> b -> b) -> b -> Tree a -> c
foldTree f g z (Node x t) = f x $ foldr (g . foldTree f g z) z t
                       -- = f x . foldr g z . map (foldTree f g z) $ t

--      1
--     / \
--    2   3
--   / \   \
--  4   5   6

t = Node 1 [Node 2 [Node 4 [], Node 5 []],
            Node 3 [Node 6 []]]

result = foldTree (+) (*) 1 t   -- ((+)%(*)) 1 t        {-
       = 1 + product [2 + product[ 4 + product[], 5 + product[]],
                      3 + product[ 6 + product[]]]
       = 1 + product [2 + 5*6, 3 + 7]
       = 321                                            -}

   -- product [a,b,c,...,n] = foldr (*) 1 [a,b,c,...,n]
   --     = a * (b * (c * ... (n * 1) ... ))

```

This is the `redtree` ("*red*uce *tree*") function from John Hughes paper, "Why Functional Programming Matters", where it is presented in a more verbose formulation.
> 7 votes

# Answer
# F#, 70 61

```
let rec r f g a (Node(n,c))=f n<|List.foldBack(g<<r f g a)c a

```

Not going to win the contest, but I think you can't get less with F#
> 5 votes

# Answer
# Haskell, 35

```
data T a=a:>[T a]
h(%)g a(e:>l)=e%foldr(g.h(%)g a)a l

```

the data type declaration is not counted.
> 4 votes

# Answer
# Prolog, 85

Logic programming

```
b([],A,A). b([E|N],A,R):-b(N,A,T),c(E,A,Y),g(Y,T,R). c((C,S),A,R):-b(S,A,T),f(C,T,R).

```

Ungolfed

```
b([],A,A).
b([Element|Siblings],A,R):-
      b(Siblings,A,RSiblings),c(Element,A,RElement),
      g(RElement,RSiblings,R).

c((Parent,Children),A,R):-
      b(Children,A,RChildren),f(Parent,RChildren,R).

```

Functions f and g

```
f(A,B,R):-R is A+B.
g(A,B,R):-R is A*B.

```

Example

```
c((1,[(2,[(4,[]),(5,[])]),(3,[(6,[])])]),1,R).

```

Try it online

Edit: correction of error signaled by Will Ness, thanks for the feedback
> 3 votes

# Answer
# Mathematica, 48

```
Head@#~#2~Fold[##3,#0[x,##2]~Table~{x,List@@#}]&

```

Example:

```
In[1]:= foldTree = Head@#~#2~Fold[##3,#0[x,##2]~Table~{x,List@@#}]&;

In[2]:= foldTree[1[2[4[], 5[]], 3[6[]]], Plus, Times, 1]

Out[2]= 321

```
> 1 votes

# Answer
# JavaScript (ES6), 62

Javascript seems functional enough, function objects can be easily defined and passed as parameters.

Define a node as an object {v:value, c:list of children}, the example tree is:

```
t={v:1, c:[{v:2,c:[{v:4},{v:5}]},{v:3,c:[{v:6}]}]}

```

(here the `c` is optional. Making `c` a mandatory field with value \[\] if there are no children, the fold function can be 5 chars shorter)

and the fold function is:

```
F=(t,a,f,g)=>{
  var r=a
  t.c && t.c.forEach(t=>r=g(r,F(t,a,f,g)))
  return f(t.v,r)
}

```

**Golfed** (there is little to golf)

```
F=(t,a,f,g,r=a)=>(t.c&&t.c.map(t=>r=g(r,F(t,a,f,g))),f(t.v,r))

```

**Test** In Firefox

```
t={v:1, c:[{v:2,c:[{v:4},{v:5}]},{v:3,c:[{v:6}]}]}

F=(t,a,f,g,r=a)=>(t.c&&t.c.map(t=>r=g(r,F(t,a,f,g))),f(t.v,r))

document.write(F(t, 1, (a,b)=>a+b, (a,b)=>a*b))
```
> 1 votes

---
Tags: code-golf, functional-programming, haskell
---


title: Gode-Colf: Wipping Flords around

text:
# Question
Title: Gode-Colf: Wipping Flords around
Here is another short string related code-golf. You must write a full program that reads a string (from STDIN or as command line arguments, as long as you specify which) and wips the flords around (flips the words around). A valid flord wip must

* Flip all of the beginning consonants. (All of the leading consonants up to the first vowel) For example, "Flipping Words" does not become "Wlipping Fords" because this only flips the first letter of each. Instead, it becomes "Wipping Flords"
* Never flip a vowel. For example, "Eating Breakfast" becomes "Breating Eakfast"
* Capitalize the first letter of each output word, and no more. For example 'Stack Overflow' becomes 'Ack Stoverflow', not 'ack StOverflow'.

Your input:

* Will have the first letter of every word capitalized, and all the other letters lowercase.
* Will consist entirely of upper and lowercase letters and spaces. There will not be any extra or leading spaces.
* There will be no strange consonant only words. For the purposes of this challenge, "y" is a consonant.
* Will have at least 2 words, each separated with a single space. If there are an even number of words, you must flip the first word with the second, the third with the fourth, the fifth with the sixth and so on and so forth. If there are an odd number of words, the same rules apply, and the last word is skipped.

Here are some sample inputs and outputs.

> Bunny Rabbit ---\> Runny Babbit
>
> Programming Puzzles And Code Golf ---\> Pogramming Pruzzles Cand Ode Golf
>
> The Cake Is A Lie ---\> Ce Thake Is A Lie ("Is" and "A" both start with a vowel, so nothing changes. "Lie" is the last odd-word-out so nothing changes)
>
> May The Best Golfer Win ---\> Thay Me Gest Bolfer Win

# Answer
# CJam, ~~62~~ 53 bytes

This is too long for now. Will add explanations once I am done golfing.

```
qS%2/{{el__"aeiou"&0=#_2$<@@>}%(a+2/Wf%~}%{s(eu\+}%S*

```

Try it online here
> 1 votes

---
Tags: code-golf, string
---