Podziękowania dla hobby Calvina za popchnięcie mojego pomysłu na wyzwanie we właściwym kierunku.

Rozważ zestaw punktów w płaszczyźnie, które nazwiemy witrynami , i skojarzmy kolor z każdą witryną. Teraz możesz pomalować całą płaszczyznę, kolorując każdy punkt kolorem najbliższego miejsca. Nazywa się to mapą Voronoi (lub diagramem Voronoi ). Zasadniczo mapy Voronoi można zdefiniować dla dowolnej metryki odległości, ale użyjemy zwykłej odległości euklidesowej r = √(x² + y²). ( Uwaga: niekoniecznie musisz wiedzieć, jak obliczać i renderować jeden z nich, aby konkurować w tym wyzwaniu).

Oto przykład ze 100 stronami:

Jeśli spojrzysz na dowolną komórkę, wszystkie punkty w tej komórce są bliżej odpowiedniego serwisu niż do dowolnego innego serwisu.

Twoim zadaniem jest przybliżenie danego obrazu za pomocą takiej mapy Voronoi. Dostaniemy obraz w dowolnym, wygodnym formacie grafiki rastrowej, jak również całkowitą N . Następnie powinieneś stworzyć do N stron i kolor dla każdej strony, tak aby mapa Voronoi oparta na tych stronach była jak najbardziej zbliżona do obrazu wejściowego.

Możesz użyć fragmentu stosu u dołu tego wyzwania, aby wyrenderować mapę Voronoi na podstawie wyników lub możesz ją renderować samodzielnie, jeśli wolisz.

Państwo może wykorzystać funkcje wbudowane lub innych firm w celu obliczenia mapy Voronoi ze zbioru stron (jeśli trzeba).

To konkurs popularności, więc wygrywa odpowiedź z największą liczbą głosów netto. Zachęca się wyborców do oceniania odpowiedzi według

• jak dobrze oryginalne obrazy i ich kolory są przybliżone.
• jak dobrze algorytm działa na różnych rodzajach obrazów.
• jak również algorytm działa dla małych N .
• czy algorytm adaptacyjnie grupuje punkty w obszarach obrazu, które wymagają więcej szczegółów.

Testuj obrazy

Oto kilka zdjęć, na których można przetestować algorytm (niektórzy z naszych zwykłych podejrzanych, niektórzy nowi). Kliknij zdjęcia, aby zobaczyć większe wersje.

Plaża w pierwszym rzędzie została narysowana przez Olivię Bell i uwzględniona za jej zgodą.

Jeśli chcesz dodatkowego wyzwania, wypróbuj Yoshi na białym tle i popraw jego linię brzucha.

Wszystkie te obrazy testowe można znaleźć w galerii imgur, z której można pobrać wszystkie w postaci pliku zip. Album zawiera również losowy schemat Voronoi jako kolejny test. Dla porównania, oto dane, które je wygenerowały .

Proszę dołączyć przykładowe diagramy dla różnych obrazów i N , np. 100, 300, 1000, 3000 (a także przypisy do niektórych odpowiednich specyfikacji komórek). Możesz używać lub pomijać czarne krawędzie między komórkami według własnego uznania (może to wyglądać lepiej na niektórych obrazach niż na innych). Nie dołączaj jednak witryn (z wyjątkiem osobnego przykładu, być może jeśli chcesz wyjaśnić, jak działa umieszczanie witryny).

Jeśli chcesz pokazać dużą liczbę wyników, możesz utworzyć galerię na imgur.com , aby zachować rozsądny rozmiar odpowiedzi. Ewentualnie umieść miniatury w swoim poście i umieść w nich linki do większych zdjęć, tak jak zrobiłem to w mojej referencji . Małe miniatury można uzyskać, dołączając sdo nazwy pliku w linku imgur.com (np. I3XrT.png-> I3XrTs.png). Ponadto, jeśli znajdziesz coś fajnego, skorzystaj z innych zdjęć testowych.

Renderer

Wklej dane wyjściowe do następującego fragmentu kodu, aby wyświetlić wyniki. Dokładny format listy nie ma znaczenia, o ile każda komórka jest określona przez 5 liczb zmiennoprzecinkowych w kolejności x y r g b, gdzie xi ysą współrzędnymi miejsca komórki, i r g bsą kanałami koloru czerwonego, zielonego i niebieskiego w zakresie 0 ≤ r, g, b ≤ 1.

Fragment zapewnia opcje określania szerokości linii krawędzi komórki i określania, czy witryny komórek powinny być pokazywane (te ostatnie głównie w celu debugowania). Pamiętaj jednak, że dane wyjściowe są ponownie renderowane tylko wtedy, gdy zmieniają się specyfikacje komórek - więc jeśli zmodyfikujesz niektóre inne opcje, dodaj spację do komórek lub czegoś innego.

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>

Ogromne podziękowania dla Raymonda Hilla za napisanie tej naprawdę miłej biblioteki JS Voronoi .

Powiązane wyzwania

• Maluj według numerów
• Photomosaics lub: Ilu programistów zajmuje wymiana żarówki?

Python + scipy + scikit-image , ważone próbkowanie płyt Poissona

Moje rozwiązanie jest dość złożone. Robię wstępne przetwarzanie obrazu, aby usunąć szum i uzyskać mapowanie, jak „interesujący” jest każdy punkt (używając kombinacji lokalnej entropii i wykrywania krawędzi):

Następnie wybieram punkty próbkowania za pomocą próbkowania dysku Poissona ze skrętem: odległość koła zależy od ciężaru, który ustaliliśmy wcześniej.

Następnie, gdy mam punkty próbkowania, dzielę obraz na segmenty voronoi i przypisuję średnią L * a * b * wartości kolorów w każdym segmencie do każdego segmentu.

Mam dużo heurystyki, a także muszę trochę matematyki, aby upewnić się, że liczba próbnych punktów jest zbliżona N. Dostaję się Ndokładnie przez nieznaczne przeregulowanie , a następnie upuszczenie niektórych punktów heurystycznie.

Pod względem czasu działania ten filtr nie jest tani , ale wykonanie poniższego zdjęcia nie zajęło więcej niż 5 sekund.

Bez zbędnych ceregieli:

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!")

Zdjęcia

Odpowiednio Nwynosi 100, 300, 1000 i 3000:

C ++

Moje podejście jest dość powolne, ale jestem bardzo zadowolony z jakości wyników, które daje, szczególnie jeśli chodzi o zachowanie krawędzi. Na przykład, oto Yoshi i Cornell Box z zaledwie 1000 stron każda:

Są dwie główne części, które sprawiają, że tyka. Pierwszy, wcielony w evaluate()funkcję, przyjmuje zestaw lokalizacji kandydujących na strony, ustawia na nich optymalne kolory i zwraca wynik dla PSNR renderowanej teselacji Voronoi względem obrazu docelowego. Kolory dla każdej witryny są określane przez uśrednienie docelowych pikseli obrazu pokrytych komórką wokół witryny. Korzystam z algorytmu Welforda, aby pomóc obliczyć zarówno najlepszy kolor dla każdej komórki, jak i wynikowy PSNR, używając tylko jednego przejścia obrazu, wykorzystując związek między wariancją, MSE i PSNR. Zmniejsza to problem do znalezienia najlepszego zestawu lokalizacji witryny bez szczególnego uwzględnienia koloru.

Druga część, wcielona w main(), próbuje znaleźć ten zestaw. Zaczyna się od losowego wyboru zestawu punktów. Następnie na każdym kroku usuwa jeden punkt (przechodząc do rundy) i testuje zestaw losowych punktów kandydujących, aby go zastąpić. Ten, który daje najwyższy PSNR w grupie, jest akceptowany i przechowywany. W efekcie strona przeskakuje do nowej lokalizacji i ogólnie poprawia obraz krok po kroku. Należy pamiętać, że algorytm celowo nie zachowuje oryginalnej lokalizacji jako kandydata. Czasami oznacza to, że skok obniża ogólną jakość obrazu. Zezwolenie na to pomaga uniknąć utknięcia w lokalnych maksimach. Daje także kryteria zatrzymania; program kończy się po przejściu określonej liczby kroków bez ulepszania najlepszego zestawu dotychczas znalezionych witryn.

Należy pamiętać, że ta implementacja jest dość podstawowa i może z łatwością zająć wiele godzin rdzenia procesora, zwłaszcza gdy liczba witryn rośnie. Ponownie oblicza pełną mapę Voronoi dla każdego kandydata, a brutalna siła testuje odległość do wszystkich miejsc dla każdego piksela. Ponieważ każda operacja wymaga usunięcia jednego punktu na raz i dodania kolejnego, rzeczywiste zmiany w obrazie na każdym kroku będą miały charakter lokalny. Istnieją algorytmy efektywnego przyrostowego aktualizowania diagramu Voronoi i wierzę, że znacznie przyspieszyłyby ten algorytm. W tym konkursie postanowiłem jednak zachować prostotę i brutalność.

Kod

#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;
 }

Bieganie

Program jest samodzielny i nie ma zewnętrznych zależności poza standardową biblioteką, ale wymaga obrazów w formacie binarnym PPM . Korzystam z ImageMagick do konwertowania obrazów do PPM, chociaż GIMP i wiele innych programów też to potrafi.

Aby go skompilować, zapisz program jako, voronoi.cppa następnie uruchom:

g++ -std=c++11 -fopenmp -O3 -o voronoi voronoi.cpp

Spodziewam się, że prawdopodobnie będzie działać w systemie Windows z najnowszymi wersjami programu Visual Studio, chociaż nie próbowałem tego. Musisz się upewnić, że kompilujesz się w C ++ 11 lub lepszej wersji i jeśli to zrobisz, masz włączony OpenMP. OpenMP nie jest absolutnie konieczne, ale bardzo pomaga w zwiększeniu tolerancji czasów wykonania.

Aby go uruchomić, zrób coś takiego:

./voronoi cornell.ppm 1000 cornell-1000.txt

Późniejszy plik będzie aktualizowany wraz z danymi witryny. Pierwszy wiersz będzie miał szerokość i wysokość obrazu, a następnie wiersze wartości x, y, r, g, b odpowiednie do kopiowania i wklejania do mechanizmu renderującego Javascript w opisie problemu.

Trzy stałe w górnej części programu pozwalają dostroić go do prędkości względem jakości. decimationCzynnikiem coarsens obraz docelowy przy ocenie zestaw stron dla koloru i PSNR. Im wyższy, tym szybciej program będzie działał. Ustawienie wartości 1 powoduje użycie obrazu w pełnej rozdzielczości. Na candidatesstałe kontrole ilu kandydatów do testowania na każdym kroku. Wyższe daje większą szansę na znalezienie odpowiedniego miejsca, do którego można skoczyć, ale powoduje spowolnienie programu. Na koniec, terminationile kroków może przejść program, nie poprawiając swoich wyników przed zakończeniem pracy. Zwiększenie go może dać lepsze wyniki, ale może nieco dłużej.

Zdjęcia

N = 100, 300, 1000 i 3000:

IDL, udoskonalenie adaptacyjne

Ta metoda jest inspirowana przez Adaptive Mesh Refinement z symulacji astronomicznych, a także powierzchnię podrejonową . Jest to zadanie, z którego szczyci się IDL, o czym będziesz wiedział dzięki dużej liczbie wbudowanych funkcji, z których mogłem korzystać. :RE

Mam kilka produktów pośrednich dla obrazu testowego Yoshi na czarnym tle, z n = 1000.

Najpierw wykonujemy skalę szarości luminancji na obrazie (używając ct_luminance) i stosujemy filtr Prewitt ( prewittpatrz wikipedia ) dla dobrego wykrywania krawędzi:

Potem następuje prawdziwa pomruk: dzielimy obraz na 4 i mierzymy wariancję w każdej ćwiartce przefiltrowanego obrazu. Nasza wariancja jest ważona rozmiarem podziału (który w tym pierwszym kroku jest równy), dzięki czemu regiony „ostre” o dużej wariancji nie stają się coraz mniejsze. Następnie używamy ważonej wariancji, aby kierować poddziały z większą liczbą szczegółów, i iteracyjnie dzielimy każdą sekcję zawierającą szczegółowe informacje na 4 dodatkowe, aż osiągniemy docelową liczbę witryn (każda podgrupa zawiera dokładnie jedną witrynę). Ponieważ dodajemy 3 witryny za każdym razem, gdy wykonujemy iterację, powstają n - 2 <= N <= nwitryny.

Zrobiłem .webm procesu podziału dla tego obrazu, którego nie mogę osadzić, ale jest tutaj ; kolor w każdej podsekcji zależy od ważonej wariancji. (Zrobiłem ten sam rodzaj wideo dla yoshi na białym tle, z odwróconą tablicą kolorów, więc staje się biały zamiast czarnego; jest tutaj .) Końcowy produkt podziału wygląda następująco:

Gdy mamy już naszą listę poddziałów, przechodzimy przez każdą z poddziałów. Ostateczną lokalizacją witryny jest pozycja minimum obrazu Prewitt, tj. Najmniej „ostry” piksel, a kolor przekroju to kolor tego piksela; oto oryginalny obraz z zaznaczonymi witrynami:

Następnie używamy wbudowanego triangulatedo obliczania triangulacji Delaunaya miejsc, a wbudowanego voronoido definiowania wierzchołków każdego wielokąta Voronoi, zanim narysujemy każdy wielokąt w buforze obrazu w odpowiednim kolorze. Na koniec zapisujemy migawkę bufora obrazu.

Kod:

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

Zadzwoń do tego przez voro_map, n, image, output_filename. Dołączyłem również wrapperprocedurę, która przeszła przez każdy obraz testowy i działała z 100, 500 i 1000 stron.

Dane wyjściowe zebrane tutaj , a oto kilka miniatur:

n = 100

n = 500

n = 1000

Python 3 + PIL + SciPy, Fuzzy k-znaczy

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)

Algorytm

Podstawową ideą jest to, że k-oznacza grupowanie naturalnie dzieli obraz na komórki Voronoi, ponieważ punkty są powiązane z najbliższym środkiem ciężkości. Musimy jednak jakoś dodać kolory jako ograniczenie.

Najpierw konwertujemy każdy piksel do przestrzeni kolorów Lab , aby lepiej manipulować kolorami.

Następnie robimy coś w rodzaju „wykrycia krawędzi biedaka”. Dla każdego piksela patrzymy na jego ortogonalnych i ukośnych sąsiadów i obliczamy średnią kwadratową różnicę koloru. Następnie sortujemy wszystkie piksele według tej różnicy, z pikselami najbardziej podobnymi do ich sąsiadów na początku listy, a pikselami najbardziej niepodobnymi do sąsiadów z tyłu (tj. Bardziej prawdopodobne, że są punktem krawędzi). Oto przykład planety, gdzie im jaśniejszy jest piksel, tym bardziej różni się od swoich sąsiadów:

(Na renderowanym wyjściu powyżej jest wyraźny wzór podobny do siatki. Według @randomra jest to prawdopodobnie spowodowane stratnym kodowaniem JPG lub kompresowaniem obrazów przez Imgur).

Następnie używamy tej kolejności pikseli do próbkowania dużej liczby punktów do zgrupowania. Używamy rozkładu wykładniczego, nadając priorytet punktom bardziej zbliżonym do krawędzi i „interesującym”.

W przypadku grupowania najpierw wybieramy Ncentroidy, losowo wybierane przy użyciu tego samego rozkładu wykładniczego, jak powyżej. Przeprowadzana jest początkowa iteracja, a dla każdej z powstałych klastrów przypisujemy średni kolor i próg wariancji kolorów. Następnie dla szeregu iteracji:

• Zbuduj triangulację centroidów w Delaunay, abyśmy mogli łatwo przesyłać zapytania sąsiadom do centroidów.
• Użyj triangulacji, aby usunąć centroidy, które są zbliżone kolorem do większości (> 4/5) swoich sąsiadów i sąsiadów razem wziętych. Wszelkie powiązane punkty próbne są również usuwane, a nowe centroidy zastępcze i punkty próbne są dodawane. Ten krok próbuje zmusić algorytm do umieszczenia większej liczby klastrów tam, gdzie potrzebne są szczegóły.
• Skonstruuj drzewo KD dla nowych centroidów, abyśmy mogli z łatwością wyszukiwać najbliższe centroidy do dowolnego punktu próbki.
• Użyj drzewa, aby przypisać każdy punkt próbki do jednego z 6 najbliższych centroidów (6 wybranych empirycznie). Środek ciężkości zaakceptuje punkt próbny tylko wtedy, gdy jego kolor mieści się w progu wariancji koloru środka ciężkości. Próbujemy przypisać każdy punkt próbki do pierwszego akceptującego centroidu, ale jeśli nie jest to możliwe, po prostu przypisujemy go do najbliższego centroidu. „Nieostrość” algorytmu pochodzi z tego kroku, ponieważ klastry mogą się nakładać.
• Przelicz centroidy.

(Kliknij, aby zobaczyć pełny rozmiar)

Na koniec próbkujemy dużą liczbę punktów, tym razem stosując równomierny rozkład. Używając innego drzewa KD, przypisujemy każdy punkt do jego najbliższego środka ciężkości, tworząc klastry. Następnie przybliżamy medianę koloru każdej gromady za pomocą algorytmu wspinania się na wzgórze, podając nasze ostateczne kolory komórek (pomysł na ten krok dzięki @PhiNotPi i @ MartinBüttner).

Notatki

Oprócz wyprowadzania pliku tekstowego dla snippet ( OUTFILE), jeśli DEBUGjest ustawiony Truena program, program również wypisze i nadpisuje powyższe obrazy. Algorytm zajmuje kilka minut dla każdego obrazu, więc jest to dobry sposób na sprawdzenie postępu, który nie wydłuża strasznie czasu pracy.

Przykładowe wyniki

N = 32:

N = 100:

N = 1000:

N = 3000:

Mathematica, Random Cells

To podstawowe rozwiązanie, więc masz pojęcie o minimum, o które cię proszę. Biorąc pod uwagę nazwę pliku (lokalnie lub jako URL) filei N do n, poniższy kod po prostu wybiera N losowych pikseli i używa kolorów znalezionych w tych pikselach. To jest naprawdę naiwne i nie działa niewiarygodnie dobrze, ale chcę, żebyście to jednak pokonali. :)

data = ImageData@Import@file;
dims = Dimensions[data][[1 ;; 2]]
{Reverse@#, data[[##]][[1 ;; 3]] & @@ Floor[1 + #]} &[dims #] & /@ 
 RandomReal[1, {n, 2}]

Oto wszystkie obrazy testowe dla N = 100 (wszystkie obrazy prowadzą do większych wersji):

Jak widać, są one zasadniczo bezużyteczne. Choć mogą mieć wartość artystyczną, w ekspresjonistyczny sposób oryginalne obrazy są ledwo rozpoznawalne.

W przypadku N = 500 sytuacja nieco się poprawiła, ale nadal istnieją bardzo dziwne artefakty, obrazy wyglądają na wyprane, a wiele komórek jest marnowanych na regiony bez szczegółów:

Twoja kolej!

Matematyka

Wszyscy wiemy, że Martin kocha Mathematicę, więc pozwól mi spróbować z Mathematica.

Mój algorytm wykorzystuje losowe punkty z krawędzi obrazu, aby utworzyć początkowy diagram voronoi. Schemat jest następnie utrwalany przez iteracyjne dopasowanie siatki za pomocą prostego filtra średniego. Daje to obrazy o wysokiej gęstości komórek w pobliżu regionów o wysokim kontraście i przyjemnych wizualnie komórek bez szalonych kątów.

Poniższe obrazy pokazują przykład procesu w akcji. Zabawa jest nieco zepsuta przez Matematicas złe antyaliasing, ale dostajemy grafikę wektorową, która musi być coś warta.

Algorytm ten, bez losowego próbkowania, można znaleźć w VoronoiMeshdokumentacji tutaj .

Obrazy testowe (100,300,1000,3000)

Kod

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

Python + SciPy + emcee

Algorytm, którego użyłem, jest następujący:

1. Zmień rozmiar zdjęć do niewielkiego rozmiaru (~ 150 pikseli)
2. Wykonaj nieostre maskowanie obrazu maksymalnych wartości kanału (pomaga to nie wychwycić zbyt mocno białych obszarów).
3. Weź wartość bezwzględną.
4. Wybierz losowe punkty z prawdopodobieństwem proporcjonalnym do tego obrazu. To wybiera punkty po obu stronach nieciągłości.
5. Udoskonal wybrane punkty, aby obniżyć funkcję kosztów. Ta funkcja to maksimum sumy kwadratowych odchyleń w kanałach (ponownie pomagając w uzyskaniu jednolitych kolorów, a nie tylko jednolitej bieli). Niewłaściwie wykorzystałem Markov Chain Monte Carlo z modułem emcee (wysoce zalecane) jako optymalizator. Procedura zostaje odrzucona, gdy po iteracjach łańcucha N nie znaleziono żadnej nowej poprawy.

Algorytm wydaje się działać bardzo dobrze. Niestety może rozsądnie działać tylko na małych obrazach. Nie miałem czasu, aby wziąć punkty Voronoi i zastosować je do większych zdjęć. W tym momencie można je udoskonalić. Mógłbym również uruchomić MCMC dłużej, aby uzyskać lepsze minima. Słabym punktem algorytmu jest to, że jest raczej drogi. Nie miałem czasu, aby przekroczyć 1000 punktów, a kilka z 1000 punktów wciąż jest udoskonalanych.

(kliknij prawym przyciskiem myszy i zobacz obraz, aby uzyskać większą wersję)

Linki do większych wersji to http://imgur.com/a/2IXDT#9 (100 punktów), http://imgur.com/a/bBQ7q (300 punktów) i http://imgur.com/a/rr8wJ (1000 punktów)

#!/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()

Nieostre zamaskowane obrazy wyglądają następująco. Punkty losowe są wybierane z obrazu, jeśli liczba losowa jest mniejsza niż wartość obrazu (znormalizowana do 1):

Opublikuję większe zdjęcia i punkty Voronoi, jeśli będę miał więcej czasu.

Edycja: Jeśli zwiększysz liczbę spacerowiczów do 100 * npts, zmień funkcję kosztu na niektóre z kwadratów odchyleń we wszystkich kanałach i poczekaj długo (zwiększając liczbę iteracji, aby wyjść z pętli do 200), możliwe jest wykonanie dobrych zdjęć za pomocą zaledwie 100 punktów:

Wykorzystanie energii obrazu jako mapy masy punktowej

W moim podejściu do tego wyzwania chciałem sposobu odwzorowania „istotności” określonego obszaru obrazu na prawdopodobieństwo, że określony punkt zostanie wybrany jako centroid Voronoi. Jednak nadal chciałem zachować artystyczny styl mozaiki Voronoi, losowo wybierając punkty obrazu. Dodatkowo chciałem operować na dużych obrazach, więc nie tracę niczego w procesie próbkowania w dół. Mój algorytm jest mniej więcej taki:

1. Dla każdego obrazu utwórz mapę ostrości. Mapa ostrości jest określona przez znormalizowaną energię obrazu (lub kwadrat sygnału wysokiej częstotliwości obrazu). Przykład wygląda następująco:

2. Wygeneruj liczbę punktów z obrazu, biorąc 70 procent z punktów na mapie ostrości i 30 procent z wszystkich innych punktów. Oznacza to, że punkty są próbkowane bardziej gęsto z fragmentów obrazu o wysokiej szczegółowości.
3. Kolor!

Wyniki

N = 100, 500, 1000, 3000