Tag Archive | splines

Splines extraction on Kopf-Lischinski algorithm, part 2


This is a sequel of a series of posts. The last one was “part 0” and this order is a bit insane, but nothing to worry about. Check the first one to understand what all this is about. I hope this will be the final post in this series.

Once upon a time… a polygon with holes appeared

The polygon-union described in the previous posts works correctly and there are no better answers than the correct answer, but these polygons were meant to be B-Splines and when I described them, I already knew that the mentioned algorithm would be insufficient.

Consider the following image.

Normal polygon (left) and bordered polygon (right)

Normal polygon (left) and bordered polygon (right)

The left polygon is the polygon that you see before and/or after the polygon-union. They are visually indistinguishable. The right polygon has an outline to allow you understand its internal representation. You can understand how the polygon-union generates a polygon with a “hole” seeing the right polygon.

If the polygon-union don’t generate visually identifiable pixels, then there shouldn’t exist any problem, but when we convert the polygon to a B-Spline, some areas of the image won’t be filled. The polygons won’t fit correctly, like show in the image below.

After B-Spline conversion

After B-Spline conversion

The solution is to create some holes. With the representation of polygons used in libdepixelize, this task is very simple and key operation to accomplish the goal is to encounter the “edges that the polygon share with itself”. I’ll explain this task further in the following paragraphs.

A holed polygon with labelled points represented with a single path

The above image has 2 pairs of hidden points, pair <5, 6>, that is equal to pair <11, 10>, and pair <12, 13>, that is equal to pair <18, 17>. Well, to extract internal polygons that represent the holes, you just iterate over the points and, for each point, try to find another point that is equal to it, then get the internal sequence and use it to construct the internal polygon.

  • Remark #1: libdepixelize will do these steps in backward order, to avoid too much moves of elements in the internal vector.
  • Remark #2: The polygon is represented in clockwise order, but the holes will be represented in counter-clockwise order, but there is no reason to normalize. You can safely ignore this (de)normalization feature.


I forgot to mention that edges not always are two points only and you should use the same “largest common edge” from the algorithm of previous posts.


Things became more complicated than I predicted and now involve some recursive functions and I needed to extend (extend != change) the algorithm. To see the whole extension, check out the Splines::_fill_holes function in libdepixelize source code and where this function is used.


Things became even more complicated than the complicated things that I didn’t anticipate. Complicated code to fix the position of T-junction nodes created a pattern where the polygon itself shared a point with itself and this pattern propagated wrongly in the “holes-detection” algorithm and afterwards.

The algorithm “3rd ed.” check if the common edge has “lenth” 1 to solve the problem.

And then, the polygon were meant to be a B-Spline

It’s a series about splines extraction and it’s fair to end with this step.

The representation of the data in three different moments

It’s evolution, baby

The algorithm is very simple. The points you get through the polygon-union algorithm in the previous posts will be the control points of quadratic Bézier curves and the interpolating/ending points of each Bézier curve will be the middle points between two adjacent control points. This is the special case where all points are smooth. The image below can help you get the idea.

Locations of the B-Spline points

Locations of the B-Spline points

If you have non-smooth points, then all you need to do is interpret them as interpolating points, instead of control points. You’ll have three interpolating points that are part of two straight lines. There is a bit more of theory to understand why the generated B-Spline is smooth, but we’re ALMOST done! In next post (a small one, I hope), I’ll detail how to get correctly positioned points for T-junctions.

Splines extraction on Kopf-Lischinski algorithm, part 1


If you don’t follow this blog and don’t know what is the Kopf-Lischinski algorithm, this small section is for you. Kopf-Lischinski algorithm is a vectorization algorithm specialized on pixel art images giving excellent results (see their supplementary material page). This “summer”, I’m working on bringing this algorithm to Inkscape in the form of a library. This post is like a “progress update status”.

Splines extraction

Well, one of the phases in the algorithm is splines extraction. This phase creates a pipeline (output from one step is the input of the next) like the following:

  1. It takes the connectivity graph
  2. It generates a generalized Voronoi diagram
  3. It groups the voronoi cells to identify visible edges
  4. It generates the splines

In a future post, I’ll explaing a fast method (already implemented in libdepixelize) to compute the Voronoi diagram exploring special properties of the input graph. In this post, I’ll explain a fast method to add polygons together exploring special properties of the generated Voronoi diagram and the data layout used in libdepixelize.

Polygon union

Polygons can be represented by vertice points. libdepixelize store them in clockwise order, with the first point being part of the polygon’s northwest/top-left and all Voronoi cells generated in step 2 are convex polygons. Another important feature of the generated Voronoi diagram is that all “connected” cells share a common edge.

The algorithm can be described in 4 major steps:

  1. Find the largest common edge
  2. Remove the in-between points of the common edge
  3. Choose one polygon and refer to it as P (choose the smallest polygon for better performance)
    1. Shift P such as P’s head and tail are points of the common edge
  4. Choose the other and refer to it as Q
    1. Remove one of the common edge’s points in Q
    2. Replace the remaining point that is part of the common edge in Q by P’s points

The Voronoi cells are iterated one by one, line by line, from the left-most cell to the right-most cell. We check if (1) we already have a group of the current color, (2) the existing group has a common edge with current Voronoi cell and, then, (3) we add the current Voronoi cell to the existing polygon. Beware that the Voronoi cells from the input are always convex polygons, but the existing polygon (that groups the Voronoi cells) might be concave.

Let’s see an example of the algorithm. In the example, we are iterating in the second line and we found an existing grouping polygon with a matching color (the two entities are connected), then we must add them togheter. The image below represents the situation we are interested in:


Polygon A is represented by the list of points [1, 2, 3, 4, 5] (common edge’s points in underline). Polygon B is represented by the list of points [6, 7, 8, 9, 10] (common edge’s points in underline). Points 4 and 8 are equal and points 5 and 7 are equal. Polygon A is the grouping polygon while polygon B is the current Voronoi cell. We shift B’s list to [8, 9, 10, 6, 7] and use it to get the final polygon [1, 2, 3, 8, 9, 10, 6, 7].

Let’s do it one more time:


Polygon A is the grouping polygon, represented by [1, 2, 3, 8, 9, 10, 6, 7]. Polygon C is the current Voronoi cell and is represented by [11, 12, 13, 14]. This time the largest common edge have 3 points and point 8/11 is the in-between point, then we remove it. The resulting A polygon is [1, 2, 3, 9, 10, 6, 7] and the resulting C polygon is [12, 13, 14]. When we replace the common edge interval in A by C, we get the final polygon, [1, 2, 12, 13, 14, 10, 6, 7]. The final polygon is show in the image below:


One last note that might be useful in later steps is the access pattern used by the algorithm. When we add a voronoi cell to the grouping polygon, there is a hot area and a cold area. The cold area is the area where there will never be a common edge. These areas always are concetrated in the same places, like exemplified by the following image:


Splines generation

Do the previous step look simple? I want to keep things simple, but there may be additional info we may want to store for each node. This (polygon-union) is the last step where we still can gather locality info without executing a deep lookup.

Valence-3 node (in red)

Valence-3 node (in red)

Let’s refer to nodes where at most two different grouping polygons share a point as valence-2 nodes. There are valence-2 nodes, valence-3 nodes and valence-4 nodes. Valence-2 nodes are always smooth and valence-4 nodes never are smooth, but valence-3 nodes vary.


Non-smooth node

Most of the points are shared by nodes of different polygons and when we have three valence-3 nodes, exactly only one of them will be smooth. We apply Kopf-Lischinski algorithm heuristics to determine which one will be and store this info for later usage. I want to play more with the code to determine the best phase to compute this data, but I’m thinking about merging Voronoi generation and node type computing in one.

Smooth point

Smooth node

The complicated part about the splines extraction on Kopf-Lischinski algorithm is the overlapping between these last steps. I wasn’t able to get the code right before “put all the cards in the table”. Now I’ll code using this guide as a reference and I’ll provide the remaining algorithm in the next post.

A bit of performance

So, Kopf-Lischinski algorithm resembles a compiler. You have several data representation types and each type offers different operations. You explore the operations each type offers and convert it to another representation until you reach the final result. In several of the type representations used by the Kopf-Lischinski algorithm, you have a matrix and you access each element and its neighbours.

A table of 3 rows and 3 columns. The elements of the first row are 1, 2 and 3. The elements of the second row are 4 violet, 5 and 6 green violet. The elements of the last row is 7, 8 and 9.

Access pattern representation

The particular implementation for libdepixelize stores and access the elements linearly like [1, 2, 3, 4, 5, 6, 7, 8, 9]. It could make good use of processor cache, but ain’t that easy. Suppose we are iterating on element 5, then we need to access all its neighbours, but only neighbours 4 and 6 may be in cache, especially in large images. This is the first problem in cache usage of the implementation, but we cannot remove this usage pattern, because it’s part of the algorithm and there is a data dependency among the elements.

What we can is reduce it, because this same pattern is used again and again over most of the passes. Another day I was reading about Xiph.org’s Daala codec and an interesting idea idea hit me: overlapping transforms. We can make each element store information about its neighbours that may not be in cache (1, 2, 3, 7, 8 and 9). This task would require more CPU time, but this cost would be mostly needed during the intialization of the first data representation only and it would be alleviated in the next passes. This change may increase memory requirements also, but the CPU would thank us for a major increase in cache hits.

Another problem with the current data access pattern is that each element store a complex object that may point to other regions of memory and add a level of indirection that can contribute to cache miss. One idea that can increase the cache hit is the one behind the Small String Optimization. This change would highly increase data locality and fits well in Kopf-Lischinski algorithm, because the complex objects stored by each element in every phase tends to have a small maximum number of subelements.

I need to combine the previous ideas and measure the impact. I have two computers and I can request access to the cluster in my university, but it’d be great to have a large database of pixel-art images ready to test. Another interesting project to use during the measurements is Stabilizer.

All that said, these are only optimizations exploring knowledge of computer architecture and I want to improve performance exploring better algorithms also. There are already a few of these, but Nathan Hurst pointed me a hint to research about Eigenvalues and I need to reserve some time to do it.

Oh, and I LOVE to work on projects where I can toy with my knowledge and projects where I can learn more and it’s so much fun to work in this project that I cannot share the joy I have right now with you guys, but I can tell this: This is the last section on this blog post, but is the first section I wrote and I only started to write the other sections after I finished this one.

Inkscape, curvas de Bézier, SVG e Python

Nesse texto, pretendo introduzir o leitor ao desenvolvimento de programas que façam uso de curvas de Bézier e, como metodologia de ensino experimental, pretendo usar conhecimentos relacionados (a ferramenta Inkscape e o padrão SVG) como auxílio.

Edição vetorial

Chegamos em um ponto onde os computadores transformaram nossos hábitos diários e estão presentes nas mais variadas partes de nossas vidas. Uma das áreas que foram modernizadas pelos computadores é a área de tratamento de imagens.

A forma mais comum utilizada para representar imagens em computadores é através de matrizes onde cada ponto contém um nível de intensidade das cores vermelho, verde e azul. Os monitores funcionam de forma parecida. O problema com essa abordagem é que se precisarmos modificar a imagem, mudando o tamanho, a posição ou o ângulo da imagem, a imagem passará por uma edição destrutiva, onde há a perda de qualidade na imagem original. Começando pelo exemplo mais simples, veja as imagens abaixo:

Small girl

Small girl

Small girl

Small girl

Note como a imagem, ao aplicarmos um “zoom”, torna-se quadriculada. É possível remover esse aspecto aplicando alguns filtros, mas não importa o filtro que escolhermos, ele irá provocar alguma perda de informação. Filtros especializados são populares, pois costumam conseguir resultados melhores para casos específicos (fotos, pixel art, …). A esse tipo de imagem, utilizamos o termo raster.

A solução que encontramos para o problema de qualidade perfeita independente de resolução foram imagens vetoriais. Imagens vetoriais são imagens que, no lugar de ter uma representação dos pontos da imagem, possuem um conjunto de descrições de curvas, que são utilizadas para gerar a imagem raster em qualquer resolução, mantendo as características da curva. A imagem abaixo ilustra bem a diferença entre imagens raster e imagens vetoriais:

Exemplo da diferença entre imagens raster e imagens vetoriais

Esse tipo de imagem ainda não é dominante, pois é mais difícil de trabalhar com elas para certos tipos de trabalho. Entretanto, para alguns tipos de trabalho (como logotipos), elas são um requisito.


Inkscape é uma ferramenta open source para trabalharmos com imagens vetoriais que vem sendo desenvolvida há bastante tempo.

O Inkscape

O Inkscape

A principal ferramenta que usamos para trabalhar com imagens vetoriais é a curva de Bézier. No Inkscape, você pode pressionar Shift + F6 para selecioná-la.

Inkscape - Curvas de Bézier

Inkscape – Curvas de Bézier

A curva de Bézier (do ponto de vista do Inkscape) é definida por nós de controle e você pode mudar certos atributos de cada nó para fazer a curva se tornar mais aberta, fechada ou cúspide.

Inkscape - edição de nós

Inkscape – edição de nós

O formato que o Inkscape usa para salvar o seu trabalho é o SVG, um padrão aberto da W3C (a mesma entidade que define vários padrões utilizados pela web) baseado em XML, que, em breve, será abordado novamente.

O foco desse artigo não é transformá-lo em um grande designer, então isso é tudo que será comentado sobre o Inkscape.


Python é uma linguagem de programação razoavelmente popular com um incrível poder expressivo que inspira vários programadores. Nesse texto ela será utilizada para a demonstração dos algoritmos e os códigos costumam ser bem legíveis, sem que nenhum treinamento especial seja necessário para o entendimento. Entretanto, utilizo alguns poucos conceitos que talvez sejam incomuns, então resolvi documentá-los aqui, para facilitar na futura compreensão de tais códigos.


Em Python, o for atua da forma que, em outras linguagens, atua o for-each. Uma função comumente utilizada quando você precisa manipular elementos através de índices é a função range, que recebe um argumento inteiro N e retorna uma lista contendo os números no intervalo [0;N). A saída de range(10), por exemplo, é [0, 1, 2, 3, 4, 5, 6, 7, 8, 9].


Em Python, existe um operador para fatiar listas (e similares). Ele possui a mesma sintaxe do operador de índice. O código a seguir resume bem algumas das operações possíveis com o mesmo:

Listas e tuplas

Em Python, os dois tipos principais para se trabalhar com cadeias são listas e tuplas. Colchetes são utilizados para criar listas, enquanto parênteses são utilizados para criar tuplas. A diferença entre listas e tuplas é que tuplas são imutáveis. Para converter entre um tipo e outro, basta passar o objeto como argumento do construtor do outro tipo.

Sobrecarregando o operador de chamada de função

Em Python, tudo é um objeto (incluindo as funções) e é possível sobrecarregar até o ponto. Para sobrecarregar o operador de chamada de função, basta definir a função-membro __call__ e, daí em diante, você poderá usar o objeto com o operador sobrecarregado como se fosse uma função.

Considerações finais sobre Python

Nesse texto, o código-fonte apresentado utiliza tuplas para representar pontos e, apesar de não haver verificações de segurança, assumimos que um ponto pode ser de qualquer dimensão (na verdade não). Isso implica que sempre utilizamos repetição para preencher os valores dos pontos.

Curvas de Bézier

Há vários modos que podemos utilizar para descrever curvas e cada abordagem reflete uma definição diferente e possui diferentes limitações. Se usarmos funções matemáticas fundamentais, como x^2 ou 2*x para representar curvas, não poderemos representar curvas que passam pela mesma linha em y mais de uma vez, uma limitação inaceitável. Utilizar funções complicadas também não é uma solução viável, pois a instabilidade numérica dos algoritmos utilizados pode estragar o resultado final. Felizmente, foram criadas as curvas de Bézier, que resolvem todos esses problemas.

Função batman

Função batman

A nossa jornada para explorar as curvas de Bézier começa com uma fórmula matemática bem simples, acompanhada de seu algoritmo em Python (por que usar pseudo-algoritmos se temos uma linguagem mais expressiva, fácil e inambígua?):


A função l(t) define uma linha que passa pelos pontos p1 e p2. No código-fonte, foi criada uma classe que permite definir uma família dessas funções, com diferentes pontos p1 e p2. Para obter um ponto presente na linha, você passa um valor t qualquer. É fácil perceber que o valor de l(0) é igual a p1 e o valor de l(1) é igual a p2.

“Curva” definida por dois pontos

Usando dois pontos e a fórmula anterior é possível definir uma linha reta. Podemos adicionar um terceiro ponto e a “mesma fórmula” para definir a curva que vai da linha (a fórmula anterior usava pontos) l1 até a linha l2. Com 3 pontos, a fórmula é transformada no seguinte:


Toda a ideia de curvas de Bézier é uma generalização desse princípio: Dada uma lista de pontos, nós queremos descrever uma curva que viaja do ponto p0 até o ponto pn guiada pelos pontos intermediários. Uma curva de Bézier composta por N pontos de controle possui grau N – 1. O código anterior se traduz na seguinte animação:

Você pode (e deveria) brincar com curvas de Bézier nesse site (que não necessita do Flash Player ou do Java).

Algoritmo de De Casteljau

Agora que sabemos o que é uma curva de Bézier, devemos começar a nos preocupar como podemos renderizá-las e esse é justamente o propósito do algoritmo de De Casteljau, um algoritmo recursivo lento, porém numericamente estável.

A propriedade explorada pelo algoritmo de De Casteljau para renderizar linhas de Bézier é que qualquer curva de Bézier pode ser dividida em duas, conectadas fim-a-fim, produzindo o mesmo resultado (lembre-se que é um algoritmo recursivo).

No código anterior, a função is_flat informa quando uma curva está indistinguível de uma reta (uma forma de heurística) e a função subdivide divide uma curva de Bézier em duas. Esse pequeno código simples é suficiente para entender as diferentes partes do algoritmo. Há uma condição para parar a recursão e duas sub-tarefas (seguindo o princípio de dividir para conquistar).

O problema será quebrado, primeiramente, com a solução para o problema de dividir uma curva de Bézier cúbica em duas. Ao tentar encontrar o ponto onde t=1/2 em uma curva de Bézier, encontramos o ponto que usaremos para “cortá-la” em duas. Na imagem a seguir, esse é o ponto azul ciano:


Existe uma afirmação de que a curva formada pelos pontos P0, P1, P2 e P3 é igual a junção da curva formada pelos pontos P0, m0, q0, B(1/2) com a curva formada pelos pontos B(1/2), q1, m2 e P3. A prova dessa afirmação existe, porém não será feita abaixo, pois ela é de pouco interesse (nem tanto…) para o nosso objetivo. Note que os pontos m0, q0, m1, q1, m2 são os pontos médios dos segmentos de retas anteriores e são facilmente computáveis. O código Python a seguir é capaz de computar esses pontos:

Note que o algoritmo serve para curvas de Bézier cúbicas (formadas por 4 pontos). Agora que temos um algoritmo para dividir as curvas de Bézier, a implementação restante para sermos capazes de utilizar o algoritmo de De Casteljau é a “heurística” que informa se a curva já está “reta” o suficiente.

Uma heurística que podemos utilizar e possui boa estabilidade numérica computa a distância, num instante t, entre onde a curva está e onde ela deveria estar caso fosse uma linha reta. Chamemos essa função de d(t), então o valor de t que queremos utilizar para a função heurística é o valor onde d(t) é maior (maximizar). Após cancelar, reduzir e otimizar várias operações, chegamos a seguinte função (que só funciona para pontos bidimensionais):


A abordagem usada até o momento requer que compartilhemos nossas imagens vetoriais através de longos textos descritivos ou códigos feito em Python, uma abordagem inviável caso queiramos compartilhar e colaborar nossos trabalhos. Para resolver esses problemas, o padrão SVG foi criado.

Como comentado anteriormente, SVG é o padrão aberto mais difundido para se trabalhar com imagens vetoriais. O SVG é baseado em XML (sendo assim, legível para humanos) que contém descrições sobre curvas de Bézier (e outros conceitos que o tornam mais expressivo e poderoso).

O elemento raiz de um SVG é o svg e seu esqueleto básico pode ser visto abaixo:

O exemplo anterior possui um comentário que indica onde colocamos o conteúdo dos nossos trabalhos, em SVG. O elemento que utilizamos para criar linhas de Bézier é o elemento path, que recebe um atributo d, capaz de controlar, de forma análoga, os movimentos que você pode fazer com um lápis. O exemplo abaixo será útil para explicar características importantes:

O documento acima possui dois novos elementos do tipo path. Eles desenham duas linhas diagonais para ajudar você a entender o sistema de coordenadas do SVG e possuem traços de preenchimento com largura e cor diferentes, para diferenciá-los mais facilmente.

Sistema de coordenadas

Sistema de coordenadas

O sistema de coordenadas é parecido com o desenho torto exemplificado acima. Você tem o “centro” no canto superior esquerdo e anotamos as coordenadas na forma x y. É diferente da forma que comumente fazemos em matemática, mas se você é um programador já deve estar acostumado com esse sistema.

Outra característica interessante para se notar é que os documentos são desenhados na ordem que foram definidos. Nesse caso, a “linha azul” irá sobrepor a “linha preta”, pois foi definida por último.

O atributo stroke define a cor de contorno (que até agora é o próprio caminho) e ela pode ser especificada de forma similar a HTML/CSS. Uma alternativa para a cor azul, usada anteriormente como blue, seria “#0000ff“. O atributo stroke-width determina a largura do contorno. A largura de contorno padrão é 1.

Para especificarmos o caminho, usamos uma linguagem que trabalha com uma sequência de instruções. A instrução M recebe uma coordenada como argumento e serve para “movimentar o lápis”. A instrução L também recebe uma coordenada como argumento e serve para fazer uma linha. Assim, essa “máquina” mantém um estado e a ordem em que você poe as instruções influencia no resultado.

Caso você tente fazer duas linhas não colineares com o que aprendeu, vai perceber que um triângulo preenchido será desenhado, pois o padrão é que a cor de preenchimento do objeto seja a cor preta. Para desativar o preenchimento, use o atributo fill com o argumento none (fill=”none”).

Hora de aprender a desenhar curvas de Bézier, o momento pelo qual estávamos esperando! O comando Q faz curvas de Bézier de quadráticas (que precisam de 3 pontos, sendo o primeiro o ponto no qual o “lápis” está), enquanto o comando C faz curvas de Bézier cúbicas (que precisam de 4 pontos). Veja o exemplo abaixo (e não fique só olhando os códigos de exemplo, teste-os):

Se você não sabia que vírgulas podem ser utilizadas como separadores de coordenadas, deve ter percebido com esse último exemplo. Outra coisa legal é que se você, em vez de digitar um comando, simplesmente colocar uma coordenada, a “máquina” entenderá como um comando L, por padrão.


Esse “trabalho” foi baseado no excelente texto de Jeremy Kun publicado no blog Math ∩ Programming, que é licenciado sob a Creative Commons Attribution-NonCommercial 3.0 Unported. Outra fonte importante foi o SVG Primer, hospedado na W3C.

Imagens vindas da Wikipédia e de outras fontes estão “corretamente” linkadas para a fonte original. Demais imagens são autoria própria (e por isso estão mais “feias”).

%d blogueiros gostam disto: