Projeto Virtual Worlds

Computação Gráfica 3D em C

Arquivos Mensais: maio 2014

Árvores de Particionamento Espacial Binarias (Binary Spatial Partitioning Trees – BSPs)

Às vezes a coisa é simples e eu faço uma confusão danada procurando significados escondidos que não estão lá… Esse foi o meu problema com as BSPs, que existem para solucionar um único problema:

Como podemos saber quais polígonos (objetos) devem ser desenhados antes de outros sem fazer muitas comparações?

O problema pode ser respondido organizando os polígonos antes de desenhá-los. Poderíamos ter uma lista com todos os polígonos e colocá-los em ordem com algum algoritmo tipo quick sort ou bubble sort, mas isso tomaria um tempo enorme. Essa aproximação ao problema torna impossível executar uma animação gráfica em tempo real, já que esses métodos levam tempo de, pelo menos, da ordem de O(N) ou O(N2), ou pior. Qual é a outra estrutura que podemos usar onde os itens são automaticamente colocados de forma ordenada e, na pesquisa, o tempo é proporcional a O(log2N)? Isso mesmo: Árvores Binárias!

Mas não podemos usar a tradicional árvore binária onde há apenas uma comparação de valores entre os nós. Nossa árvore precisa ter dois tipos de nós: Os nós contendo informações sobre o “mundo” e os nós contendo informações sobre os polígonos que serão desenhados em cada pedaço desse mundo. Precisamos dividir nosso mundo em pedaços, chamados partições (daí o nome “Particionamento Espacial”). Nossa árvore tem, então, nós (nodes) contendo planos e nós “folha” (leaf nodes), contendo os polígonos que pertencem a uma partição do mundo.

Como acredito que já tenha mostrado, planos são artefatos matemáticos onde todos os pontos contidos nele obedecem uma equação simples:

\displaystyle \vec{n}\cdot\dot{p}+d=0

Onde \vec{n} é o vetor unitário normal (perpendicular) ao plano, \dot{p} é a coordenada de um ponto na superfície do plano e d é a distância deste a origem do sistema de coordenadas à qualquer ponto do plano (por isso d é calculado pela projeção sobre o vetor normal e tem seu sinal invertido!). A operação “∙” é um produto escalar, o que denota a projeção do ponto \dot{p} sobre o vetor normal \vec{n}. O ponto \dot{p} é a variável da equação.

Você também encontra a equação de um plano, em livros de geometria, escrito dessa maneira:

\displaystyle Ax+By+Cz+D=0

Onde A, B e C são os componentes X, Y e Z do vetor unitário normal \vec{n}, as variáveis x, y e z são as coordenadas do ponto \dot{p} e D é a distância da projeção de \dot{p} sobre \vec{n}, em relação à origem… Com isso, podemos determinar se um outro ponto qualquer \dot{q} está na frente, sobre ou atrás do plano definido pela equação. E, só para constar, prefiro a equação vetorial… Acho naus intuitiva e elegante!

Voltando à BSP, se cada nó que não é leaf contém a equação do plano que particiona o mundo (ou o sub-mundo, digamos assim) em dois, podemos determinar se um ponto \dot{q} qualquer está na frente ou atrás do plano em questão. Basta colocar as coordenadas do ponto \dot{q} no lugar do ponto \dot{p} e verificar o valor final da equação do plano… Se o valor for zero, o ponto encontra-se no plano, se for maior que zero, está na frente, caso contrário, atrás. Ao determinar isso podemos decidir se vamos continuar percorrendo a árvore do lado da frente ou de trás do plano, repetindo o teste com os outros planos que encontrarmos até acharmos um nó leaf.

O que queremos? Desenhar as regiões em ordem de proximidade…

A primeira coisa que faremos é determinar o nó leaf onde o nosso ponto \dot{p} está contido (agora uso o ponto \dot{p} como o ponto da câmera, não o da equação)… Considere esse simples mapa, onde as partes pretas são paredes (que considerarei como tendo espessura muito fina, para facilitar a explicação) e as áreas delimitadas pelos planos (linhas vermelhas) e paredes fronteiriças são “setores” definidos pelas partições (o parte pontilhada do “plano” c tem uma explicação que explicarei logo):
map
O plano a divide o mapa na região A e nas regiões B, C e D. Ele divide o mapa em duas partes…

O plano b divide o lado da frente do plano a nas regiões B, de um lado, e C e D, do outro.
Por sua vez, o plano c divide C e D, que eram uma região só, em duas partes.

Se escolhermos o plano b para ser o nó raiz de nossa árvore, teremos:

tree
Onde, o que está à esquerda de um nó é tudo aquilo que está atrás dele e, por consequência, o que está à direita, está na frente. Assim, o plano a está atrás do plano b e região D está na frente do plano b (mas, atrás do plano c). Eu escolhi o plano b como raiz por que sabia, de antemão, que obteria uma árvore perfeitamente balanceada. A escolha do plano divisor raiz é essencial para obter uma árvore “ótima”. Mas isso pode ser conseguido usando-se um algoritmo de árvore balanceada…

Agora que temos nossa BSP montada, tudo o que temos que fazer é percorrê-la em pre-order (raiz primeiro, depois os nós filhos). A cada visita a um nó não leaf (os nós pretos) temos que fazer o teste para saber se nossa câmera está na frente ou atrás do plano do nó. Suponha que nossa câmera esteja na parte inferior da região A.

Estamos no nó raiz: A câmera encontra-se em frente ou atrás do plano b? Para saber isso basta usar a equação do plano, lá em cima… No caso, a câmera está atrás do plano b, então vamos, recursivamente, para o nó da esquerda e fazemos a pergunta de novo, mas para o plano a. Ora, a câmera está atrás do plano a, então vamos para a esquerda, de novo, e encontramos o nó leaf contendo a região A. Este nó contém todos os polígonos que estão na região, então os desenhamos todos, exceto aqueles que falham nos testes de visibilidade (frunstum culling e aqueles cujos vetores normais estejam na mesma direção do vetor da câmera, por exemplo). Já que o nó A não tem filhos, a rotina recursiva retornará, fazendo com que o nó B, filho do nó do plano a, seja visitado. Ao retornar, o nó de c será visitado e um novo teste com o plano é feito. Nossa câmera pode estar ou não atrás do plano c (e eis o motivo do tracejado no mapa). Já nossa câmera está na parte inferior da região A, ela estará atrás do plano c, fazendo com que o nó da região D seja visitada e depois o nó da região C.

Note que as regiões visitadas seguiram a ordem A, B, D e C.

Se a câmera estivesse na região superior de B, a ordem seria: B, A, C e D (faça o teste!). Se estivesse na região D, teríamos: D, C, B e A… Repare que a ordem das regiões segue um padrão útil: São regiões ordenadas de acordo com a proximidade de onde câmera está!

E é só para isso que serve o particionamento espacial: Colocar as regiões em ordem! A partir daí, podemos usar outros algoritmos (que usam a BSP) para determinar quais superfícies escondidas podem ser removidas (Hidden Surface Removal) e qual conjunto de polígonos são potencialmente visíveis (Potentially Visible Set), mas isso eu deixo para outra ocasião… O código para percorrer a árvore BSP, uma vez que ela tenha sido montada, é mais ou menos esse:

struct bspnode_s {
  struct bspnode_s *front, *back;
  float plane[4]; // equação do plano, para nós não "folha".
  struct objectlist_s objlist; // Lista de objetos na região (nó folha).
};

void bsp_traverse(struct bspnode_s *node)
{
  struct bspnode_s *f, *b;

  if (node->front == NULL && node->back == NULL)
    drawobjlist(node->objlist);
  else
  {
    // Assume que a câmera esteja na frente do plano.
    f = node->front;
    b = node->back;

    // A câmera está atrás do plano? Então trocamos a ordem de pesquisa!
    if (PlaneGetDistance(node->plane, camera) < 0.0f)
    {
      f = b;
      b = node->front;
    }

    // Recursivamente visita os nós filhos, na ordem especificada no teste acima.
    bspnode_traverse(f);
    bspnode_traverse(b);
  }
}

O código que “compila” a BSP eu também deixo para outra ocasião. Ele é meio complicado porque envolve alguns detalhes interessantes, como aquele do plano pontilhado e outros, como, determinação de conjunto convexo de polígonos, sem contar com o balanceamento da árvore. Sobre esse último item, em breve escreverei um artigo para o Bit Is Myth sobre Red-Black Trees.

O motivo de ter descartado o formato Wavefront

Estou trabalhando na carga de objetos construídos a partir do Blender. Ao invés de usar o formato binário do editor preferi usar um formato mais “amigável” para importar os objetos e, talvez, criar meu próprio formato binário. Um dos formatos mais simples de lidar é o da WaveFront Technologies (a mesma empresa que criou o Maya)… Só que vi que o exporter default do Blender para esse formato não me dá aquilo que desejo:

  1. Aparentemente os agrupamentos de objetos não são exportados;
  2. exporter não otimiza os meshes para triangle strips ou triangle fans;

Um exemplo é este: Criei dois cubos (e “triangulei” as faces!). O da direita tem o nome de Cube1 e o da esquerda, Cube2. Coloquei ambos os objetos num grupo chamado Group1 e exportei para .OBJ. Abaixo você vê um grab da tela do Blender e o arquivo cubes.obj gerado:

Cubes

Eis o arquivo cubes.obj:

# Blender v2.63 (sub 0) OBJ File: 'untitled.blend'
# www.blender.org
mtllib cubes.mtl
o Cube2
v -1.000000 -1.000000 -1.000000
v -1.000000 -1.000000 1.000000
v -3.000000 -1.000000 1.000000
v -3.000000 -1.000000 -1.000000
v -1.000000 1.000000 -0.999999
v -1.000001 1.000000 1.000001
v -3.000000 1.000000 1.000000
v -3.000000 1.000000 -1.000000
vn 0.000000 -1.000000 0.000000
vn -0.000000 1.000000 0.000000
vn 1.000000 -0.000000 0.000000
vn -0.000000 -0.000000 1.000000
vn -1.000000 -0.000000 -0.000000
vn 0.000000 0.000000 -1.000000
vn 1.000000 0.000000 0.000001
usemtl Material
s off
f 1//1 2//1 3//1
f 5//2 8//2 6//2
f 1//3 5//3 2//3
f 2//4 6//4 3//4
f 3//5 7//5 8//5
f 5//6 1//6 4//6
f 4//1 1//1 3//1
f 8//2 7//2 6//2
f 5//7 6//7 2//7
f 6//4 7//4 3//4
f 4//5 3//5 8//5
f 8//6 5//6 4//6
o Cube1
v 3.000000 -1.000000 -1.000000
v 3.000000 -1.000000 1.000000
v 1.000000 -1.000000 1.000000
v 1.000000 -1.000000 -1.000000
v 3.000000 1.000000 -0.999999
v 2.999999 1.000000 1.000001
v 1.000000 1.000000 1.000000
v 1.000000 1.000000 -1.000000
usemtl Material
s off
f 9//1 10//1 11//1
f 13//2 16//2 14//2
f 9//3 13//3 10//3
f 10//4 14//4 11//4
f 11//5 15//5 16//5
f 13//6 9//6 12//6
f 12//1 9//1 11//1
f 16//2 15//2 14//2
f 13//7 14//7 10//7
f 14//4 15//4 11//4
f 12//5 11//5 16//5
f 16//6 13//6 12//6

Algumas coisas que você pode observar:

  1. O grupo Group1 não foi exportado;
  2. A lista de vértices e vetores normais (e se tivesse, a lista de coordenadas de texturas) são compartilhadas pelos objetos − e isso é bom! − Mas, não há quaisquer informações sobre strips ou fans;

Eu ainda não sei se formatos mais “modernos” como COLLADA e X3D fazem esses tipos de otimizações. Mas, estou estudando ambos para ver se abandono de vez o formato Wavefront.

Abaixo, temos o código (parcial) de como montar um objeto lendo o arquivo .OBJ. Deixei o tratamento de texturas e outros “comandos” fora desde código para ilustrar a facilidade de lidar com o formato… Começe olhando para LoadOBJ().

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <setjmp.h>
#include <malloc.h>

/* Tamanho máximo de uma linha do arquivo .OBJ. */
#define MAX_LINE_SIZE 255

/* Mantém a lista de vértices e normais. */
static float *vertices = NULL;
static float *normals = NULL;
static int numvertices = 0;
static int numnormals = 0;

/* Estrutura que contém os índices dos triangulos. */
struct triangle_s {
  float vertice[3];
  float normal[3];
};

/* Estrutura de um "objeto". */
struct object_s {
  char name[32];

  struct triangle_s *triangles;
  int numtriangles;
};

/* Mantém a lista de "objetos" */
struct object_s *objects = NULL;
int numobjects = 0;

/* Realoca espaço para uma "lista" */
static void *grow_array(void **pp, size_t itemsize, size_t numitems)
{
  void *tmpp;

  if ((tmpp = realloc(*pp, itemsize * numitems)) != NULL)
    *pp = tmpp;

  return tmpp;
}

/* Lê uma linha do arquivo, retirando '\n' do final e
   retornando um ponteiro para o primeiro caracter não 'vazio'. */
static char *ReadLine(FILE *f, char *buffer, size_t max_size)
{
  char *cp;

  if ((cp = fgets(buffer, max_size, f)) != NULL)
  {
    char *tmpcp;

    cp += strspn(buffer, " \t\r");
    if ((tmpcp = strchr(cp, '\n')) != NULL)
      *tmpcp = '';
  }
  return cp;
}

static void Vector3Copy(float *dest, const float *src);
static int GetAllVerticesAndNormals(FILE *fin);
static int GetObjectsAndTriangles(FILE *fin);

/* Carrega .OBJ, retorna -1 em caso de erro ou 0 se ok. */
int LoadOBJ(char *filename)
{
  FILE *fin;
  jmp_buf jb; /* Usando setjmp() para tratar erro. */

  if ((fin = fopen(filename, "rt")) == NULL)
    return -1;

  if (setjmp(jb))
  {
    int i;

    fclose(fin);

    /* Joga todos as listas "no lixo" e retorna -1. */
    free(vertices); vertices = NULL; numvertices = 0;
    free(normals); normals = NULL; numnormals = 0;
    for (i = 0; i < numobjects; i++)
    {
      free((objects + i)->triangles);
      (objects + i)->triangles = NULL;
      (objects + i)->numtriangles = 0;
    } 
    free(objects); objects = NULL; numobjects = 0;

    return -1;
  }

  /* A leitura do arquivo .OBJ é feita em dois passos. Primeiro,
     lemos todos os vértices, vetores e vetores normais. */
  if (GetAllVerticesAndNormals(fin) < 0)
    longjmp(jb, 1);

  /* Passo 2: Lemos os objetos e suas faces, associando-as com os
     objetos. */
  if (GetObjectsAndTriangles(fin) < 0)
    longjmp(jb, 1);

  /* Neste ponto não precisamos mais dos vetores vertices e normals. */
  free(vertices); vertices = NULL; numvertices = 0;
  free(normals); normals = NULL; numnormals = 0;

  fclose(fin);

  /* Agora temos a lista de objetos em "objects" pronta para ser
     colocada em vertex buffers e enviadas ao OpenGL via DrawArrays(). */
  return 0;
}

/* Lê todos os vértices e vetores normais do arquivo. */
static int GetAllVerticesAndNormals(FILE *fin)
{
  char line[MAX_LINE_SIZE+1];
  char *cp;

  fseek(fin, 0L, SEEK_SET);

  while ((cp = ReadLine(fin, line, MAX_LINE_SIZE)) != NULL)
  {
    char *cp2;

    /* Isola o comando dos seus parâmetros. */
    if ((cp2 = strpbrk(cp, " \t")) != NULL)
      *cp2++ = '';

    /* Temos um vértice? */
    if (strcasecmp(cp, "v") == 0)
    {
      float *vp;

      if ((vp = grow_array((void **)&vertices, 3*sizeof(float), 
                  numvertices+1)) == NULL)
        return -1;
 
      vp += 3 * numvertices++; // Aponta para último item.

      /* Eu DEVERIA testar o retorno de sscanf(), mas... What the hell,
         isto é só um exemplo, certo? */
      sscanf(cp2, "%f %f %f", vp, vp + 1, vp + 2);
    }

    /* Temos um vertor normal? */
    if (strcasecmp(cp, "vn") == 0)
    {
      float *np;

      if ((np = grow_array((void **)&normals, 3*sizeof(float), 
                  numnormals+1)) == NULL)
        return -1;
 
      np += 3 * numnormals++; // Aponta para último item.

      /* Eu DEVERIA testar o retorno de sscanf(), mas... What the hell,
         isto é só um exemplo, certo? */
      sscanf(cp2, "%f %f %f", np, np + 1, np + 2);
    }
  } 

  return 0;
}

/* Lê objetos e "faces" */
static int GetObjectsAndTriangles(FILE *fin)
{
  char line[MAX_LINE_SIZE+1];
  char *cp;
  struct object_s *op = NULL;
  int i;

  fseek(fin, 0L, SEEK_SET);

  while ((cp = ReadLine(fin, line, MAX_LINE_SIZE)) != NULL)
  {
    char *cp2;

    /* Isola o comando dos seus parâmetros. */
    if ((cp2 = strpbrk(cp, " \t")) != NULL)
      *cp2++ = '';

    /* Temos um novo objeto? */
    if (strcasecmp(cp, "o") == 0)
    {
      if ((op = grow_array((void **)&objects, sizeof(struct object_s), 
                  numobjects + 1)) == NULL)
        return -1;

      op += numobjects++;

      strcpy(op->name, cp2);
      op->triangles = NULL;
      op->numtriangles = 0;
    }

    /* Temos um triangulo? */
    if (strcasecmp(cp, "f") == 0)
    {
      int indexes[6];
      struct triangle_s *tp;

      /* Temos que ter alocado o último objeto, pelo menos! */
      if (op == NULL)
        return -1;

      if ((tp = grow_array((void **)&op->triangles, 
                           3*sizeof(struct triangle_s), 
                           op->numtriangles + 1)) == NULL)
        return -1;

      tp += 3*op->numtriangles++;

      /* Eu DEVERIA testar o retorno de sscanf(), mas... What the hell,
         isto é só um exemplo, certo? */

      /* Obtém os índices. */
      sscanf(cp2, "%d//%d %d//%d %d//%d",
        indexes, indexes + 1,
        indexes + 2, indexes + 3,
        indexes + 4, indexes + 5);

      /* Usa os índices (que são baseados em 1, não em 0!) para
         copiar os vetores. */
      for (i = 0; i < 3; i++, tp++)
      {
        Vector3Copy(tp->vertice, vertices + 3*indexes[0+3*i] - 1);
        Vector3Copy(tp->normal, normals + 3*indexes[1+3*i] - 1);
      }
    }
  }

  return 0;
}

/* Copia um vetor para outro. */
static void Vector3Copy(float *dest, const float *src)
{
  dest[0] = src[0];
  dest[1] = src[1];
  dest[2] = src[2];
}

PS: O código acima é um exemplo e não for depurado!

Voltando atrás…

Voltei atrás e retirei a maioria das otimizações usando SSE. Especialmente com relação aos vetores. O motivo é que colocar o componente w nos vetores, embora facilite o uso de SSE, consome 4 bytes adicionais para cada vetor. Num objeto com 10000 vértices isso dá cerca de 400 kB adicionais.

As rotinas de vetores agora tém nomes como Vector3Normalize() e aceitam ponteiros para float.

Manipulação da câmera (viewpoint).

Existe uma maneira simples de lidar com a analogia da câmera, em OpenGL. Na biblioteca GLU (OpenGL Utilities) existe a função gluLookAt, que toma 3 vetores:

void gluLookAt(float eyeX, float eyeY, float eyeZ,
               float centerX, float centerY, float centerZ,
               float upX, float upY, float upZ);

O problema é que os vetores eye e up tem que ser perpendiculares.

Mas, se mudarmos o ponto de visão (o vetor “eye”), o vetor “up”, com toda certeza, não será mais perpendicular ao primeiro. E como não temos o terceiro vetor perpendicular como referência, não podemos usar um produto vetorial para calcular o novo vetor “up”. Poderíamos colocar um quarto vetor na jogada (além do vetor “eye” − que dá o ponto para onde estamos olhando; o vetor “center” − que nos dá a posição no espaço da câmera; e do vetor “up” − que nos diz onde é que nossa câmera interpreta “para cima”), mas isso só iria causar mais complicações e diminuiria a performance, já que teríamos que lidar com mais um vetor. Existe um jeito mais fácil.

Dado um vetor “eye” e um vetor “up”, não perpendicular, podemos simplesmente calcular a projeção do vetor “up” num vetor unitário baseado em “eye” e subtrair do vetor “up” essa projeção, lembrando de normalizar o vetor “up” no fim das contas. “Projeção” é “produto escalar”. A coisa fica assim:

static void RecalcUpVector(float *up, const float *eye)
{
  float w[3], v[4];

  /* w = normalize(eye); */
  Vector3Normalize(w, eye);

  /* up = up - (w . up)w */
  Vector3Scale(v, w, Vector3Dot(w, up));
  Vector3Sub(up, up, v);
  Vector3Normalize2(up);
}

Todo vetor tem 3 componentes. Para facilitar o entendimento eis como funciona com vetores 2D: Todo vetor 2D tem componentes X e Y. Ao projetar o vetor “up”, não perpendicular, sobre o vetor unitário colinear ao vetor “eye” temos um vetor apenas com o tamanho X de “up”. Sabendo que \vec{v}_{up} é dado por:

\vec{v}_{up} = \vec{v}_{up_x} + \vec{v}_{up_y}

vector

Para obter o vetor \vec{v}_{up_y}, que é perpendicular ao vetor unitário \vec{w} basta subtrair o vetor \vec{v}_{up_x} (a projeção de \vec{v}_{up} sobre \vec{w}) do próprio \vec{v}_{up}.

Ok, toda essa história de “projeção” e “normalização” pode ser resumida assim: A rotina acima substrai o componente do vetor “up” da parte que se encontra “sobre” o vetor “eye” e certifica-se que o novo vetor “up” tem tamnho igual a 1.

Se nossa câmera é definida somente com esses 3 vetores (eye, center e up), então o código para ajustar a matriz de transformação do sistema de coordenadas VIEW fica:

void SetupViewspaceMatrix(struct camera_s *camera)
{
  RecalcUpVector(camera->up, camera->eye);

  gluLookAt(camera->eye[0], camera->eye[1], camera->eye[2],
            camera->center[0], camera->center[1], camera->center[2],
            camera->up[0], camera->up[1], camera->up[2]);
}

Simples assim…

Para mover a cãmera, basta atualizar o vetor “center” (que é uma coordenada). Para rotacionar a câmera podemos usar a função Vector3RotateAroundAxis() e depois recalcular o vetor “up”. É claro que temos que ter cuidado para que os vetores “eye” e “up” não sejam colineares!