Desenhar uma imagem como um mapa de Voronoi

170

Créditos ao Hobbies de Calvin por empurrar minha ideia de desafio na direção certa.

Considere um conjunto de pontos no plano, que chamaremos de sites , e associe uma cor a cada site. Agora você pode pintar o plano inteiro colorindo cada ponto com a cor do site mais próximo. Isso é chamado de mapa de Voronoi (ou diagrama de Voronoi ). Em princípio, os mapas de Voronoi podem ser definidos para qualquer métrica de distância, mas simplesmente usaremos a distância euclidiana usual r = √(x² + y²). ( Observação: você não precisa necessariamente saber como calcular e renderizar um deles para competir nesse desafio.)

Aqui está um exemplo com 100 sites:

insira a descrição da imagem aqui

Se você olhar para qualquer célula, todos os pontos nessa célula estarão mais próximos do site correspondente do que em qualquer outro site.

Sua tarefa é aproximar uma determinada imagem com um mapa de Voronoi. Você recebe a imagem em qualquer formato gráfico de varredura conveniente, além de um N inteiro . Você deve produzir até N sites e uma cor para cada site, de modo que o mapa Voronoi baseado nesses sites se assemelhe à imagem de entrada o mais próximo possível.

Você pode usar o snippet de pilha na parte inferior deste desafio para renderizar um mapa Voronoi a partir de sua saída ou pode renderizá-lo você mesmo, se preferir.

Você pode usar funções internas ou de terceiros para calcular um mapa Voronoi a partir de um conjunto de sites (se necessário).

Como é um concurso de popularidade, ganha a resposta com o maior número de votos. Os eleitores são incentivados a julgar as respostas

  • quão bem as imagens originais e suas cores são aproximadas.
  • quão bem o algoritmo funciona em diferentes tipos de imagens.
  • quão bem o algoritmo funciona para N pequeno .
  • se o algoritmo agrupa adaptativamente pontos nas regiões da imagem que exigem mais detalhes.

Imagens de teste

Aqui estão algumas imagens para testar seu algoritmo (alguns de nossos suspeitos comuns, outros novos). Clique nas imagens para versões maiores.

Great Wave Ouriço de praia Cornell Saturno Urso marrom Yoshi Mandrill Nebulosa do Caranguejo Garoto de Geobits Cascata Grito

A praia na primeira fila foi desenhada por Olivia Bell e incluída com sua permissão.

Se você quer um desafio extra, tente Yoshi com um fundo branco e acerte a linha da barriga.

Você pode encontrar todas essas imagens de teste nesta galeria de imagens, onde é possível fazer o download de todas elas como um arquivo zip. O álbum também contém um diagrama aleatório de Voronoi como outro teste. Para referência, aqui estão os dados que os geraram .

Inclua diagramas de exemplo para uma variedade de imagens diferentes e N , por exemplo, 100, 300, 1000, 3000 (além de pastas para algumas das especificações de célula correspondentes). Você pode usar ou omitir bordas pretas entre as células como achar melhor (isso pode parecer melhor em algumas imagens do que em outras). Não inclua os sites (exceto em um exemplo separado, talvez se você quiser explicar como o posicionamento do site funciona, é claro).

Se você deseja mostrar um grande número de resultados, pode criar uma galeria no imgur.com , para manter o tamanho das respostas razoáveis. Como alternativa, coloque miniaturas em sua postagem e faça links para imagens maiores, como fiz na resposta de referência . Você pode obter as pequenas miniaturas anexando so nome do arquivo no link imgur.com (por exemplo, I3XrT.png-> I3XrTs.png). Além disso, fique à vontade para usar outras imagens de teste, se encontrar algo legal.

Renderer

Cole sua saída no seguinte snippet de pilha para renderizar seus resultados. O formato exato da lista é irrelevante, desde que cada célula seja especificada por 5 números de ponto flutuante na ordem x y r g b, onde xe ysão as coordenadas do site da célula e r g bos canais de cores vermelho, verde e azul no intervalo 0 ≤ r, g, b ≤ 1.

O trecho fornece opções para especificar uma largura de linha das bordas da célula e se os sites da célula devem ou não ser exibidos (o último principalmente para fins de depuração). Mas observe que a saída só é renderizada novamente quando as especificações da célula são alteradas - portanto, se você modificar algumas das outras opções, adicione um espaço às células ou algo assim.

Créditos maciços a Raymond Hill por escrever esta realmente agradável biblioteca JS Voronoi .

Desafios relacionados

Martin Ender
fonte
5
@frogeyedpeas Ao olhar para os votos que você recebe. ;) Este é um concurso de popularidade. Não há necessariamente a melhor maneira de fazer. A idéia é que você tente fazer o melhor possível e os votos refletirão se as pessoas concordam que você fez um bom trabalho. Há uma certa quantidade de subjetividade nelas, é certo. Dê uma olhada nos desafios relacionados que vinculei ou neste . Você verá que geralmente há uma grande variedade de abordagens, mas o sistema de votação ajuda as melhores soluções a chegarem ao topo e a decidir sobre um vencedor.
Martin Ender
3
Olivia aprova as aproximações de sua praia enviadas até agora.
Alex A.
3
@AlexA. Devon aprova algumas das aproximações de seu rosto apresentadas até agora. Ele não é um grande fã de qualquer um dos n = 100 versões;)
Geobits
1
@ Geobits: Ele entenderá quando for mais velho.
Alex A.
1
Aqui está uma página sobre uma técnica de pontilhado à base de voronoi centróide . Uma boa fonte de inspiração (a dissertação de mestrado relacionada tem uma boa discussão sobre possíveis melhorias no algoritmo).
Job

Respostas:

112

Amostra de disco de Poisson ponderada em Python + scipy + scikit-image

Minha solução é bastante complexa. Eu faço um pré-processamento na imagem para remover o ruído e obter um mapeamento de quão interessante é cada ponto (usando uma combinação de entropia local e detecção de borda):

Depois, escolho os pontos de amostragem usando a amostragem de disco de Poisson com um giro: a distância do círculo é determinada pelo peso que determinamos anteriormente.

Depois que tenho os pontos de amostragem, divido a imagem em segmentos voronoi e atribuo a média L * a * b * dos valores de cores dentro de cada segmento a cada segmento.

Tenho muitas heurísticas e também preciso fazer um pouco de matemática para garantir que o número de pontos de amostra esteja próximo N. Fico Nexatamente superando um pouco e depois removendo alguns pontos com uma heurística.

Em termos de tempo de execução, esse filtro não é barato , mas nenhuma imagem abaixo levou mais de 5 segundos para ser criada.

Sem mais delongas:

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

Imagens

Respectivamente, Né 100, 300, 1000 e 3000:

abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc
abc abc abc abc

orlp
fonte
2
Eu gosto disso; parece um pouco com vidro fumado.
BobTheAwesome
3
Eu brinquei um pouco com isso, e você obtém melhores resultados, principalmente para as imagens de triângulo baixo, se você substituir o denoise_bilatteral por um denoise_tv_bregman. Ele gera mais patches no seu denoising, o que ajuda.
LKlevin #
@LKlevin Qual peso você usou?
orlp 18/05/19
Eu usei 1.0 como o peso.
LKlevin
65

C ++

Minha abordagem é bastante lenta, mas estou muito feliz com a qualidade dos resultados que ela fornece, principalmente com relação à preservação de arestas. Por exemplo, aqui está Yoshi e a Cornell Box com apenas 1000 sites cada:

Existem duas partes principais que o fazem funcionar. O primeiro, incorporado na evaluate()função, pega um conjunto de locais candidatos, define as cores ideais e retorna uma pontuação para o PSNR do mosaico Voronoi renderizado versus a imagem de destino. As cores de cada site são determinadas pela média dos pixels da imagem de destino cobertos pela célula ao redor do site. Uso o algoritmo de Welford para ajudar a calcular a melhor cor para cada célula e o PSNR resultante, usando apenas uma única passagem sobre a imagem, explorando a relação entre variação, MSE e PSNR. Isso reduz o problema a encontrar o melhor conjunto de locais do site sem levar em consideração a cor.

A segunda parte, então, incorporada main(), tenta encontrar esse conjunto. Começa escolhendo um conjunto de pontos aleatoriamente. Em seguida, em cada etapa, ele remove um ponto (round-robin) e testa um conjunto de pontos candidatos aleatórios para substituí-lo. O que produz o PSNR mais alto do grupo é aceito e mantido. Efetivamente, isso faz com que o site salte para um novo local e geralmente melhora a imagem bit a bit. Observe que o algoritmo intencionalmente não retém o local original como candidato. Às vezes, isso significa que o salto diminui a qualidade geral da imagem. Permitir que isso aconteça ajuda a evitar ficar preso nos máximos locais. Também fornece um critério de parada; o programa termina após executar um certo número de etapas sem melhorar o melhor conjunto de sites encontrado até o momento.

Observe que essa implementação é bastante básica e pode facilmente levar horas de tempo no núcleo da CPU, especialmente à medida que o número de sites aumenta. Ele recalcula o mapa completo de Voronoi para todos os candidatos e a força bruta testa a distância de todos os sites de cada pixel. Como cada operação envolve remover um ponto de cada vez e adicionar outro, as alterações reais na imagem em cada etapa serão razoavelmente locais. Existem algoritmos para a atualização incremental eficiente de um diagrama de Voronoi e acredito que eles dariam a esse algoritmo uma tremenda aceleração. Para este concurso, no entanto, escolhi manter as coisas simples e com força bruta.

Código

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

Corrida

O programa é independente e não possui dependências externas além da biblioteca padrão, mas exige que as imagens estejam no formato binário PPM . Eu uso o ImageMagick para converter imagens em PPM, embora o GIMP e muitos outros programas também possam fazê-lo.

Para compilá-lo, salve o programa como voronoi.cppe, em seguida, execute:

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

Espero que provavelmente funcione no Windows com versões recentes do Visual Studio, embora ainda não tenha tentado. Você quer ter certeza de que está compilando com o C ++ 11 ou melhor e o OpenMP ativado, se o fizer. O OpenMP não é estritamente necessário, mas ajuda muito a tornar os tempos de execução mais toleráveis.

Para executá-lo, faça algo como:

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

O arquivo posterior será atualizado com os dados do site. A primeira linha terá a largura e a altura da imagem, seguidas pelas linhas dos valores x, y, r, g, b adequados para copiar e colar no renderizador Javascript na descrição do problema.

As três constantes na parte superior do programa permitem ajustá-lo para velocidade versus qualidade. O decimationfator aumenta a imagem de destino ao avaliar um conjunto de sites para cores e PSNR. Quanto mais alto, mais rápido o programa será executado. A configuração para 1 usa a imagem em resolução máxima. A candidatesconstante controla quantos candidatos testar em cada etapa. Maior oferece uma chance maior de encontrar um bom local para o qual pular, mas torna o programa mais lento. Finalmente, terminationé quantas etapas o programa pode executar sem melhorar sua produção antes de sair. Aumentá-lo pode dar melhores resultados, mas leva um pouco mais de tempo.

Imagens

N = 100, 300, 1000 e 3000:

Boojum
fonte
1
Isso deveria ter ganho IMO - muito melhor que o meu.
orlp 31/05
1
@orlp - Obrigado! Para ser justo, você postou o seu muito antes, e ele corre muito mais rapidamente. A velocidade conta!
Boojum
1
Bem, o meu não é realmente uma resposta do mapa voronoi :) É um algoritmo de amostragem realmente muito bom, mas transformar pontos de amostra em sites voronoi claramente não é o ideal.
orlp
55

IDL, refinamento adaptável

Este método é inspirado no Refinamento de Malha Adaptativa a partir de simulações astronômicas e também na superfície de subdivisão . Esse é o tipo de tarefa em que a IDL se orgulha, que você poderá contar pelo grande número de funções internas que pude usar. : D

Eu produzi alguns dos intermediários para a imagem de teste yoshi de fundo preto, com n = 1000.

Primeiro, executamos uma escala de cinza da luminância na imagem (usando ct_luminance) e aplicamos um filtro Prewitt ( prewittconsulte a Wikipedia ) para uma boa detecção de borda:

abc abc

Depois vem o trabalho grunhido real: subdividimos a imagem em 4 e medimos a variação em cada quadrante da imagem filtrada. Nossa variância é ponderada pelo tamanho da subdivisão (que neste primeiro passo é igual), para que regiões "nervosas" com alta variância não sejam subdivididas cada vez menores. Em seguida, usamos a variação ponderada para segmentar subdivisões com mais detalhes e subdividimos iterativamente cada seção rica em detalhes em 4 adicionais, até atingirmos o número alvo de sites (cada subdivisão contém exatamente um site). Como adicionamos três sites cada vez que iteramos, terminamos com n - 2 <= N <= nsites.

Fiz uma .webm do processo de subdivisão para esta imagem, que não consigo incorporar, mas está aqui ; a cor em cada subseção é determinada pela variação ponderada. (Fiz o mesmo tipo de vídeo para o yoshi de fundo branco, para comparação, com a tabela de cores invertida para que ela fique em branco em vez de preto; está aqui .) O produto final da subdivisão é semelhante a este:

abc

Depois de termos nossa lista de subdivisões, passamos por cada subdivisão. A localização final do site é a posição do mínimo da imagem Prewitt, ou seja, o pixel menos "nervoso", e a cor da seção é a cor desse pixel; aqui está a imagem original, com sites marcados:

abc

Em seguida, usamos o interno triangulatepara calcular a triangulação de Delaunay dos sites e o interno voronoipara definir os vértices de cada polígono de Voronoi, antes de desenhar cada polígono em um buffer de imagem em sua respectiva cor. Por fim, salvamos um instantâneo do buffer de imagem.

abc

O código:

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

Ligue para isso via voro_map, n, image, output_filename. Também incluí um wrapperprocedimento, que passou por cada imagem de teste e foi executado com 100, 500 e 1000 sites.

Saída coletada aqui e aqui estão algumas miniaturas:

n = 100

abc abc abc abc abc abc abc abc abc abc abc abc abc

n = 500

abc abc abc abc abc abc abc abc abc abc abc abc abc

n = 1000

abc abc abc abc abc abc abc abc abc abc abc abc abc

sirpercival
fonte
9
Eu realmente gosto do fato de que esta solução coloca mais pontos em áreas mais complexas, que é a intenção, e a diferencia das outras neste momento.
Alexander-brett
sim, a idéia de pontos agrupados por detalhes é o que me levou a um refinamento adaptativo #
1150 sirpercival
3
Explicação muito elegante, e as imagens são impressionantes! Eu tenho uma pergunta - parece que você obtém imagens muito diferentes quando Yoshi está em um fundo branco, onde temos algumas formas estranhas. O que pode estar causando isso?
BrainSteel 17/05
2
@BrianSteel I imaginar os contornos são apanhados como áreas de alta variância e focado em desnecessariamente, e, em seguida, outros verdadeiramente áreas de alto detalhe tem menos pontos atribuídos por causa disso.
Doppelgreener 17/05
@BrainSteel, acho que o doppel está certo - há uma forte margem entre a borda preta e o fundo branco, o que exige muitos detalhes no algoritmo. eu não tenho certeza se isso é algo que eu posso (ou, mais importante, deve ) fixar ...
sirpercival
47

Python 3 + PIL + SciPy, k difuso

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)

O algoritmo

A idéia principal é que o agrupamento de k-significa particione naturalmente a imagem nas células Voronoi, uma vez que os pontos estão ligados ao centróide mais próximo. No entanto, precisamos adicionar de alguma forma as cores como uma restrição.

Primeiro, convertemos cada pixel no espaço de cores laboratório , para melhor manipulação das cores.

Então fazemos uma espécie de "detecção de ponta do pobre homem". Para cada pixel, examinamos seus vizinhos ortogonais e diagonais e calculamos a diferença média quadrática da cor. Em seguida, classificamos todos os pixels por essa diferença, com os pixels mais semelhantes aos vizinhos na frente da lista e os pixels mais diferentes dos vizinhos na parte de trás (ou seja, é mais provável que seja um ponto de extremidade). Aqui está um exemplo para o planeta, onde quanto mais brilhante o pixel, mais diferente ele é de seus vizinhos:

insira a descrição da imagem aqui

(Existe um padrão claro de grade na saída renderizada acima. De acordo com @randomra, isso provavelmente se deve à codificação JPG com perdas ou à imgur de compactação das imagens.)

Em seguida, usamos essa ordem de pixels para provar um grande número de pontos a serem agrupados. Usamos uma distribuição exponencial, priorizando pontos que são mais parecidos com arestas e "interessantes".

insira a descrição da imagem aqui

Para o agrupamento, primeiro escolhemos Ncentróides, escolhidos aleatoriamente usando a mesma distribuição exponencial acima. Uma iteração inicial é executada e, para cada um dos clusters resultantes, atribuímos uma cor média e um limite de variação de cores. Em seguida, para várias iterações, nós:

  • Construa a triangulação de Delaunay dos centróides, para que possamos consultar facilmente os vizinhos.
  • Use a triangulação para remover os centróides com cores próximas à maioria (> 4/5) de seus vizinhos e vizinhos do vizinho combinados. Quaisquer pontos de amostra associados também são removidos e novos centróides de substituição e pontos de amostra são adicionados. Esta etapa tenta forçar o algoritmo a colocar mais clusters onde os detalhes são necessários.
  • Construa uma árvore kd para os novos centróides, para que possamos consultar facilmente os centróides mais próximos de qualquer ponto de amostra.
  • Use a árvore para atribuir cada ponto de amostra a um dos 6 centróides mais próximos (6 escolhidos empiricamente). Um centróide só aceitará um ponto de amostra se a cor do ponto estiver dentro do limite de variação de cores do centróide. Tentamos atribuir cada ponto de amostra ao primeiro centróide aceitante, mas se isso não for possível, basta atribuí-lo ao centróide mais próximo. A "imprecisão" do algoritmo vem dessa etapa, pois é possível que os clusters se sobreponham.
  • Recompute os centróides.

insira a descrição da imagem aqui

(Clique para ampliar)

Finalmente, amostramos um grande número de pontos, desta vez usando uma distribuição uniforme. Usando outra árvore kd, atribuímos cada ponto ao centróide mais próximo, formando aglomerados. Em seguida, aproximamos a cor mediana de cada cluster usando um algoritmo de escalada, fornecendo as cores finais das células (idéia para esta etapa graças a @PhiNotPi e @ MartinBüttner).

insira a descrição da imagem aqui

Notas

Além de emitir um arquivo de texto para o trecho ( OUTFILE), se DEBUGestiver definido como Trueo programa, também emitirá e sobrescreverá as imagens acima. O algoritmo leva alguns minutos para cada imagem, por isso é uma boa maneira de verificar o progresso que não adiciona muito ao tempo de execução.

Saídas de amostra

N = 32:

insira a descrição da imagem aqui

N = 100:

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

N = 1000:

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

N = 3000:

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Sp3000
fonte
1
Eu realmente gosto do quão bem seus Yoshis brancos ficaram.
Max
26

Mathematica, Células Aleatórias

Esta é a solução de linha de base, para que você tenha uma idéia do mínimo que estou pedindo a você. Dado o nome do arquivo (localmente ou como um URL) filee N in n, o código a seguir simplesmente seleciona N pixels aleatórios e usa as cores encontradas nesses pixels. Isso é realmente ingênuo e não funciona incrivelmente bem, mas quero que vocês superem isso depois de tudo. :)

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

Aqui estão todas as imagens de teste para N = 100 (todas as imagens possuem link para versões maiores):

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Como você pode ver, estes são essencialmente inúteis. Embora possam ter algum valor artístico, de uma maneira expressionista, as imagens originais são quase irreconhecíveis.

Para N = 500 , a situação melhora um pouco, mas ainda existem artefatos muito estranhos, as imagens parecem desbotadas e muitas células são desperdiçadas em regiões sem detalhes:

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Sua vez!

Martin Ender
fonte
Eu não sou um bom programador, mas meu Deus, essas imagens são lindas. Excelente ideia!
Faraz Masroor
Alguma razão para Dimensions@ImageDatae não ImageDimensions? Você pode evitar o lento ImageDatacompletamente usando PixelValue.
-campion
@ 2012rcampion Sem motivo, eu simplesmente não sabia que nenhuma das funções existia. Posso corrigir isso mais tarde e também alterar as imagens de exemplo para os valores N recomendados.
Martin Ender
23

Mathematica

Todos sabemos que Martin ama o Mathematica, então deixe-me tentar com o Mathematica.

Meu algoritmo usa pontos aleatórios das bordas da imagem para criar um diagrama voronoi inicial. O diagrama é então pré-modificado por um ajuste iterativo da malha com um filtro médio simples. Isso fornece imagens com alta densidade celular perto de regiões de alto contraste e células visualmente agradáveis ​​sem ângulos malucos.

As imagens a seguir mostram um exemplo do processo em ação. A diversão é um pouco estragada pelo Antialiasing ruim do Mathematicas, mas temos gráficos vetoriais, que devem valer alguma coisa.

Esse algoritmo, sem a amostragem aleatória, pode ser encontrado na VoronoiMeshdocumentação aqui .

insira a descrição da imagem aqui insira a descrição da imagem aqui

Imagens de teste (100,300,1000,3000)

Código

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]}]]
   ];
pata
fonte
Bom trabalho para um primeiro post! :) Você pode tentar a imagem de teste Voronoi com 32 células (esse é o número de células na própria imagem).
Martin Ender
Obrigado! Eu acho que meu algoritmo terá um ótimo desempenho neste exemplo. As sementes serão inicializados nas bordas celulares e a recursão não vai fazer isso muito melhor;)
pata
Apesar da taxa mais lenta de convergência para a imagem original, acho que seu algoritmo produz um resultado muito artístico! (como uma versão melhorada das obras de Georges Seurat). Bom trabalho!
Neizod
Você também pode obter vidrado olhando cores polígono interpolados mudando suas linhas finaisGraphics@Table[ Append[mp, VertexColors -> RGBColor /@ ImageValue[img, First[mp]]], {mp, MeshPrimitives[voronoiRelaxed, 2]}]
histogramas
13

Python + SciPy + emcee

O algoritmo que usei é o seguinte:

  1. Redimensionar imagens para um tamanho pequeno (~ 150 pixels)
  2. Faça uma imagem sem máscara dos valores máximos de canal (isso ajuda a não capturar regiões brancas com muita força).
  3. Tome o valor absoluto.
  4. Escolha pontos aleatórios com uma probabilidade proporcional a esta imagem. Isso escolhe pontos de ambos os lados das descontinuidades.
  5. Refine os pontos escolhidos para reduzir uma função de custo. A função é o máximo da soma dos desvios quadrados dos canais (novamente ajudando a polarizar cores sólidas e não apenas o branco sólido). Eu usei mal o Markov Chain Monte Carlo com o módulo emcee (altamente recomendado) como um otimizador. O procedimento falha quando nenhuma nova melhoria é encontrada após as iterações da cadeia N.

O algoritmo parece funcionar muito bem. Infelizmente, ele só pode ser executado sensivelmente em imagens pequenas. Não tive tempo de pegar os pontos Voronoi e aplicá-los às imagens maiores. Eles podem ser refinados neste momento. Eu também poderia ter executado o MCMC por mais tempo para obter melhores mínimos. O ponto fraco do algoritmo é que ele é bastante caro. Não tive tempo de aumentar além de 1000 pontos e algumas das imagens de 1000 pontos ainda estão sendo refinadas.

(clique com o botão direito do mouse e visualize a imagem para obter uma versão maior)

Miniaturas para 100, 300 e 1000 pontos

Os links para versões maiores são http://imgur.com/a/2IXDT#9 (100 pontos), http://imgur.com/a/bBQ7q (300 pontos) e http://imgur.com/a/rr8wJ (1000 pontos)

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

Imagens mascaradas sem nitidez se parecem com as seguintes. Pontos aleatórios são selecionados na imagem se um número aleatório for menor que o valor da imagem (normatizado para 1):

Imagem de Saturno mascarada sem nitidez

Vou postar imagens maiores e os pontos Voronoi, se tiver mais tempo.

Edit: Se você aumentar o número de walkers para 100 * npts, altere a função de custo para alguns quadrados dos desvios em todos os canais e aguarde um longo tempo (aumentando o número de iterações para interromper o loop para 200), é possível criar boas imagens com apenas 100 pontos:

Imagem 11, 100 pontos Imagem 2, 100 pontos Imagem 4, 100 pontos Imagem 10, 100 pontos

xioxox
fonte
3

Usando a energia da imagem como um mapa de pontos-peso

Na minha abordagem para esse desafio, eu queria uma maneira de mapear a "relevância" de uma área de imagem específica com a probabilidade de um ponto específico ser escolhido como um centróide de Voronoi. No entanto, eu ainda queria preservar a sensação artística do mosaico de Voronoi escolhendo aleatoriamente pontos de imagem. Além disso, eu queria operar imagens grandes, para não perder nada no processo de redução de amostras. Meu algoritmo é mais ou menos assim:

  1. Para cada imagem, crie um mapa de nitidez. O mapa de nitidez é definido pela energia normalizada da imagem (ou pelo quadrado do sinal de alta frequência da imagem). Um exemplo é assim:

Mapa de nitidez

  1. Gere vários pontos da imagem, obtendo 70% dos pontos no mapa de nitidez e 30% de todos os outros pontos. Isso significa que os pontos são amostrados com mais densidade de partes da imagem com alto detalhe.
  2. Cor!

Resultados

N = 100, 500, 1000, 3000

Imagem 1, N = 100 Imagem 1, N = 500 Imagem 1, N = 1000 Imagem 1, N = 3000

Imagem 2, N = 100 Imagem 2, N = 500 Imagem 2, N = 1000 Imagem 2, N = 3000

Imagem 3, N = 100 Imagem 3, N = 500 Imagem 3, N = 1000 Imagem 3, N = 3000

Imagem 4, N = 100 Imagem 4, N = 500 Imagem 4, N = 1000 Imagem 4, N = 3000

Imagem 5, N = 100 Imagem 5, N = 500 Imagem 5, N = 1000 Imagem 5, N = 3000

Imagem 6, N = 100 Imagem 6, N = 500 Imagem 6, N = 1000 Imagem 6, N = 3000

Imagem 7, N = 100 Imagem 7, N = 500 Imagem 7, N = 1000 Imagem 7, N = 3000

Imagem 8, N = 100 Imagem 8, N = 500 Imagem 8, N = 1000 Imagem 8, N = 3000

Imagem 9, N = 100 Imagem 9, N = 500 Imagem 9, N = 1000 Imagem 9, N = 3000

Imagem 10, N = 100 Imagem 10, N = 500 Imagem 10, N = 1000 Imagem 10, N = 3000

Imagem 11, N = 100 Imagem 11, N = 500 Imagem 11, N = 1000 Imagem 11, N = 3000

Imagem 12, N = 100 Imagem 12, N = 500 Imagem 12, N = 1000 Imagem 12, N = 3000

Imagem 13, N = 100 Imagem 13, N = 500 Imagem 13, N = 1000 Imagem 13, N = 3000

Imagem 14, N = 100 Imagem 14, N = 500 Imagem 14, N = 1000 Imagem 14, N = 3000

mprat
fonte
14
Você se importaria em) incluir o código-fonte usado para gerar isso eb) vincular cada miniatura à imagem em tamanho real?
Martin Ender