Skip to content

Dessiner une image comme une carte de Voronoï

Cette division a été testée par des experts afin que la véracité de ces nouvelles soit garantie.

Solution :

Python + scipy + scikit-image, échantillonnage par disques de Poisson pondérés.

Ma solution est plutôt complexe. Je fais un peu de prétraitement sur l'image pour enlever le bruit et obtenir une cartographie à quel point chaque point est 'intéressant' (en utilisant une combinaison d'entropie locale et de détection des bords) :

Ensuite, je choisis des points d'échantillonnage en utilisant l'échantillonnage par disque de Poisson avec une torsion : la distance du cercle est déterminée par le poids que nous avons déterminé précédemment.

Ensuite, une fois que j'ai les points d'échantillonnage, je divise l'image en segments de voronoï et j'attribue à chaque segment la moyenne L*a*b* des valeurs de couleur à l'intérieur de chaque segment.

J'ai beaucoup d'heuristiques, et je dois aussi faire un peu de mathématiques pour m'assurer que le nombre de points d'échantillonnage est proche de... N. J'obtiens N exactement en dépassant légèrement, puis en abandonnant quelques points avec une heuristique.

En termes de temps d'exécution, ce filtre n'est pas... bon marché, mais aucune image ci-dessous n'a pris plus de 5 secondes à réaliser.

Sans plus attendre :

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

Respectivement N est de 100, 300, 1000 et 3000 :

abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc
abcabcabcabc

C++

Mon approche est assez lente, mais je suis très satisfait de la qualité des résultats qu'elle donne, notamment en ce qui concerne la préservation des bords. Par exemple, voici Yoshi et la boîte de Cornell avec seulement 1000 sites chacun :

Il y a deux parties principales qui le font fonctionner. La première, incarnée par le evaluate() prend un ensemble d'emplacements de sites candidats, leur attribue les couleurs optimales et renvoie un score pour le PSNR de la tessellation de Voronoï rendue par rapport à l'image cible. Les couleurs pour chaque site sont déterminées en faisant la moyenne des pixels de l'image cible couverts par la cellule autour du site. J'utilise l'algorithme de Welford pour calculer à la fois la meilleure couleur pour chaque cellule et le PSNR résultant en utilisant un seul passage sur l'image en exploitant la relation entre la variance, le MSE et le PSNR. Cela réduit le problème à celui de trouver le meilleur ensemble d'emplacements de sites sans considération particulière pour la couleur.

La deuxième partie alors, concrétisée par main(), essaie de trouver cet ensemble. Elle commence par choisir un ensemble de points au hasard. Puis, à chaque étape, elle supprime un point (en faisant un tour de table) et teste un ensemble de points candidats aléatoires pour le remplacer. Celui qui donne le PSNR le plus élevé du groupe est accepté et conservé. En fait, cela fait sauter le site vers un nouvel emplacement et améliore généralement l'image bit par bit. Notez que l'algorithme fait intentionnellement pas conserver l'emplacement original comme candidat. Parfois, cela signifie que le saut diminue la qualité globale de l'image. Permettre que cela se produise permet d'éviter de rester bloqué dans des maxima locaux. Cela donne également un critère d'arrêta; le programme se termine après avoir franchi un certain nombre d'étapes sans améliorer le meilleur ensemble de sites trouvés jusqu'à présent.

Notez que cette implémentation est assez basique et peut facilement prendre des heures de temps CPU-core, en particulier lorsque le nombre de sites augmente. Elle recalcule la carte de Voronoï complète pour chaque candidat et teste par force brute la distance à tous les sites pour chaque pixel. Comme chaque opération implique la suppression d'un point à la fois et l'ajout d'un autre, les modifications réelles de l'image à chaque étape seront assez locales. Il existe des algorithmes permettant de mettre à jour efficacement et de manière incrémentielle un diagramme de Voronoï et je pense qu'ils permettraient à cet algorithme de gagner énormément de temps. Pour ce concours, cependant, j'ai choisi de garder les choses simples et la force brute.

Code

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

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 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 &sites) {
     auto counts = vector(sites.size());
     auto variance = vector(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(sites[index].x), 0), target.width - 1);
             auto y = min(max(static_cast(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 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{0.0f, 1.0f};
     auto target = read(argv[1]);
     auto sites = vector{};
     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;
 }

Exécution de

Le programme est autonome et n'a pas de dépendances externes au-delà de la bibliothèque standard, mais il nécessite que les images soient au format binaire PPM. J'utilise ImageMagick pour convertir les images en PPM, bien que GIMP et pas mal d'autres programmes puissent le faire aussi.

Pour le compiler, enregistrez le programme sous voronoi.cpp puis exécutez-le :

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

Je pense que cela fonctionnera probablement sous Windows avec des versions récentes de Visual Studio, bien que je n'ai pas essayé. Vous devrez vous assurer que vous compilez avec C++11 ou mieux et OpenMP activé si vous le faites. OpenMP n'est pas strictement nécessaire, mais il aide beaucoup à rendre les temps d'exécution plus tolérables.

Pour l'exécuter, faites quelque chose comme :

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

Le fichier ultérieur sera mis à jour avec les données du site au fur et à mesure. La première ligne aura la largeur et la hauteur de l'image, suivie de lignes de valeurs x, y, r, g, b appropriées pour être copiées et collées dans le moteur de rendu Javascript dans la description du problème.

Les trois constantes en haut du programme vous permettent de le régler pour la vitesse par rapport à la qualité. Le site decimation grossit l'image cible lors de l'évaluation d'un ensemble de sites pour la couleur et le PSNR. Plus il est élevé, plus le programme s'exécute rapidement. En le réglant sur 1, on utilise l'image à pleine résolution. Le facteur candidates constant contrôle le nombre de candidats à tester à chaque étape. Un nombre plus élevé donne une meilleure chance de trouver un bon endroit où sauter, mais rend le programme plus lent. Enfin, termination est le nombre de pas que le programme peut franchir sans améliorer sa sortie avant de s'arrêter. L'augmenter peut donner de meilleurs résultats mais le rendre marginalement plus long.

Images

N = 100, 300, 1000, et 3000 :

IDL, raffinement adaptatif

Cette méthode s'inspire de l'Adaptive Mesh Refinement des simulations astronomiques, et aussi de la surface de Subdivision. C'est le genre de tâche dont IDL s'enorgueillit, ce que vous pourrez constater par le grand nombre de fonctions intégrées que j'ai pu utiliser. 😀

J'ai sorti quelques-uns des intermédiaires pour l'image de test du yoshi à fond noir, avec n = 1000.

Tout d'abord, nous effectuons une échelle de gris de luminance sur l'image (avec la fonction ct_luminance), et nous appliquons un filtre Prewitt (prewitt, voir wikipedia) pour une bonne détection des bords :

abcabc

Puis vient le vrai travail de grogne : on subdivise l'image en 4, et on mesure la variance de chaque quadrant dans l'image filtrée. Notre variance est pondérée par la taille de la subdivision (qui, dans cette première étape, est égale), de sorte que les régions "nerveuses" avec une variance élevée ne soient pas subdivisées de plus en plus petites. Ensuite, nous utilisons la variance pondérée pour cibler les subdivisions avec plus de détails, et nous subdivisons itérativement chaque section riche en détails en 4 autres, jusqu'à ce que nous atteignions notre nombre cible de sites (chaque subdivision contient exactement un site). Puisque nous ajoutons 3 sites chaque fois que nous itérons, nous obtenons n - 2 <= N <= n sites.

J'ai fait un .webm du processus de subdivision pour cette image, que je ne peux pas intégrer, mais c'est elle...e; la couleur dans chaque sous-section est déterminée par la variance pondérée. (J'ai fait le même genre de vidéo pour le yoshi à fond blanc, à titre de comparaison, avec la table des couleurs inversée pour qu'elle aille vers le blanc au lieu du noir ; c'est ici). Le produit final de la subdivision ressemble à ceci :

abc

Une fois que nous avons notre liste de subdivisions, nous passons en revue chaque subdivision. L'emplacement final du site est la position du minimum de l'image de Prewitt, c'est-à-dire le pixel le moins "edgy", et la couleur de la section est la couleur de ce pixel ; voici l'image originale, avec les sites marqués :

abc

Ensuite, nous utilisons la fonction intégrée triangulate intégrée pour calculer la triangulation de Delaunay des sites, et la fonction intégrée voronoi pour définir les sommets de chaque polygone de Voronoï, avant de dessiner chaque polygone dans un tampon d'image dans sa couleur respective. Enfin, nous sauvegardons un instantané du tampon d'image.

abc

Le 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

Appelez ceci via voro_map, n, image, output_filename. J'ai inclus un wrapper procédure également, qui est passée par chaque image de test et a fonctionné avec 100, 500 et 1000 sites.

Sortie recueillie ici, et voici quelques vignettes :

n = 100

abcabcabcabcabcabcabcabcabcabcabcabcabc

n = 500

abcabcabcabcabcabcabcabcabcabcabcabcabc

n = 1000

abcabcabcabcabcabcabcabcabcabcabcabcabc

Section des critiques et des évaluations

Si vous aimez le sujet, vous pouvez laisser un essai sur ce que vous pensez de cette écriture.



Utilisez notre moteur de recherche

Ricerca
Generic filters

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.