Projeto Virtual Worlds

Computação Gráfica 3D em C

Arquivos da Categoria: Matemática

Mudança de base e “bump mapping”

Eis mais uma demonstração do poder das matrizes. Especialmente quanto falamos de mudança de base do sistema de coordenadas (ou, no jargão da álgebra linear: do subespaço vetorial).

Dê uma olhada nas 3 texturas abaixo:

Textura original, mapa de normais e textura modificada.

Textura original, mapa de normais e textura modificada.

A textura da esquerda é a original, parece com um muro de pedra, mas é “achatada” (se você rotacioná-la verá o efeito). A textura do meio, azulada, é, na verdade, um mapa com os vetores normais de cada pixel da textura, usada para causar o efeito final, o da textura da direita, onde o “rejunte” está  num nível abaixo do das pedras.

Cada componente da textura do meio equivale ao comprimento de um vetor numa base bem definida: R corresponde ao vetor associado ao eixo X; G ao Y; e B ao Z. A imagem é azulada porque o vetor Z é o mesmo do vetor normal original. O vetor normal original é o mais proeminente. Os outros dois são usados para modificar o vetor normal original. Por isso as bordas das pedras, nessa textura, têm cores tendendo ao violeta (componente R grande) ou ao “ciano” (componente G grande).

Visualmente, esses vetores correspondem ao diagrama abaixo:

Visualizando os vetores tangente e binormal.

Visualizando os vetores tangente e binormal.

Essa textura central é baseada, ou seja, tem como base do sistema de coordenadas, 3 vetores unitários: O vetor normal original e dois outros, chamados de tangente e binormal (ou bitangente). Tudo o que um shader tem que fazer para obter o vetor original é saber qual é a base e aplicar a mudança de base abaixo:

\displaystyle \vec{N}' = \begin{bmatrix} \vec{t}^T \\ \vec{b}^T \\ \vec{n}^T \end{bmatrix} \cdot \vec{N}

Onde \vec{t} é o vetor tangente, \vec{b} é o binormal e \vec{n} é o normal; obtidos do vertex buffer, junto com os outros atributos do vértice. \vec{N} são as 3 coordenadas obtidas da textura do meio e \vec{N}' é o vetor normal recalculado. Como quem fará esse calculo é o fragment shader, o vetor \vec{N} é obtido da textura via um Sampler. Ou seja, cada pixel da textura é um vetor. Mas, diferente da figura anterior, cada um dos vetores “tangentes” e o vetor “normal” são definidos para cada um dos vértices.

Note que a equação acima, na verdade, é uma conversão de base do sistema de coordenadas original, onde o vetor normal está para outro, baseado na superfície do objeto e \vec{N} “perturba” o vetor normal original… Ainda existe um passo intermediário: Os valores de \vec{N} precisam estar na faixa entre [-1,1], ao invés de [0,1], como acontece com os valores RGB numa textura. Isso pode ser facilmente feito escalonando \vec{N} duas vezes (multiplicando por 2) e subtraindo 1.

Calculando a iluminação:

Ainda não chegamos lá, mas um modelo para calcular iluminação (especialmente iluminação “ambiente”) é através do produto escalar entre o vetor normal e o vetor diretor da fonte de luz. Já que estamos “dando um solavanco” (bumping) os vetores normais de cada pixel, a impressão que ficará é a de que a superfície do objeto texturizado é, de fato, parecida com um muro, com rejunte e tudo. Pode-se enxergar esses “solavancos” na figura simplificada abaixo:

Antes e depois do "solavanco".

Antes e depois do “solavanco”.

O problema do bump mapping:

Nem tudo é perfeito… Bump Mapping tem um problema: A ilusão de profundidade só pode ser mantida em alguns casos. Se a superfície estiver voltada diretamente para a câmera, a ilusão ficará perfeita, mas rotacione a superfície e o encanto se quebra. A figura abaixo ilustra:

Bump Mapping e Parallax Mapping, lado a lado.

Bump Mapping e Parallax Mapping, lado a lado.

Do lado esquerdo temos um cubo texturizado e usando a técnica descrita acima. Do lado direito, temos o mesmo cubo, usando uma técnica diferente, mais complexa, chamada Steep Parallax Mapping.

Em ambos os casos temos o mesmo cubo (liso) com a mesma textura aplicada, mas com diferentes tipos de mapas de vetores normais, usando diferentes técnicas. Descreverei Steep Parallax Mapping assim que entendê-la melhor…

Isso não significa que bump mapping seja ruim… Dê uma olhada, por exemplo, no jogo Watch Dogs. Especialmente nos teclados que aparecem no jogo, sobre mesas… Eles são colocados lá por bump mapping e, ás vezes, dá para perceber o efeito acima… mas, só às vezes! :)

Explicando transformações lineares melhor do que eu poderia fazê-lo…

Veja esse playlist: Clique aqui

Uma outra interpretação de Matrix4x4LookAt

Não causa estranheza que a matriz de transformação usada por Matrix4x4LookAt seja transversa? Num de meus momentos de “meditação” sobre o assunto percebi outra possível interpretação para tal fato. Lembram-se que o produto escalar de um vetor \vec{v} por um vetor unitário \vec{u} nos dá a projeção do primeiro sobre o segundo? E que um produto escalar é definido por:

\displaystyle\vec{v}\cdot\vec{u}=v_x\cdot u_x+v_y\cdot u_y+v_z\cdot u_z

Assim, faz todo sentido obter as projeções de um vetor \vec{v} sobre os eixos do novo sistema de coordenadas através da multiplicação da matrix transversa M^T, não é? Com isso você obtem os novos valores de (x,y,z) do vetor no novo sistema de coordenadas:

\vec{{v_x}'}=\vec{v}\cdot\vec{v_{side}} \\  \vec{{v_y}'}=\vec{v}\cdot\vec{v_{up}} \\  \vec{{v_z}'}=-\vec{v}\cdot\vec{v_{fw}}

Ou, de forma matricial:

\displaystyle \vec{v'} = \begin{bmatrix} \vec{side}^T \\ \vec{up}^T \\ -\vec{fw}^T \\ 1 \end{bmatrix} \cdot \vec{v}

Simples assim.

Explicando Matrix4x4LookAt

No código de Matrix4x4LookAt, no último artigo, eu disse que ao multiplicar a matriz transposta composta pelos vetores sideup e -forward obtemos uma “reorientação” do sistema de coordenadas view. Mas, por quê?

A coisa toda pode ser entendida assim: Considere que uma transformação geométrica sobre um vetor v qualquer é dada por:

\displaystyle \vec{v'} = M\vec{v}

Onde \vec{v'} é o vetor transformado (reposicionado), M é a matriz de transformação e \vec{v} é o vetor original. Se usarmos uma matriz de identidade, temos (exemplo em 2D):

\displaystyle \vec{v'} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \cdot \vec{v}

\displaystyle \begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} x \\ y \end{bmatrix}

Uma matriz nada mais é do que uma coleção de vetores. No caso da matriz de identidade os vetores são [1 0]T e [0 1]T:

\displaystyle \vec{i} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\,\quad \vec{j} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}\,\quad \begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix} \vec{i} & \vec{j} \end{bmatrix} \cdot \begin{bmatrix} x \\ y \end{bmatrix}

O primeito vetor é o alinhado com o eixo x e o segundo, com o eixo y. A transformação acima é chamada “identidade” porque, no fim das contas, \vec{v'} será idêntico ao vetor \vec{v}.

Acontece que também podemos fazer uma transformação “contrária” (M^{-1}\cdot\vec{v'}=\vec{v}), ou seja, se a matriz estava sendo multiplicada do lado direito, ela passa pro lado esquerdo “dividindo”, mas como a divisão de matrizes ou vetores e matrizes não é definida, temos que inverter a matriz. Qualquer um que já tenha estudado algebra linear sabe que calcular uma matriz inversa é uma tarefa um pouco chata. Envolve alguns passos. Só que podemos usar um atalho em alguns casos: Quando uma matriz é composta de vetores unitários e perpendiculares entre si, então sua inversa é idêntica à sua transposta. Este é o caso de sistemas de coordenadas que nos interessam. O sitema de cordenada view é composto de 3 vetores “diretores”, perpendiculares entre si e o “novo” sistema de coordenadas – o da câmera – também (sideup-forward)…

No contexto da função Matrix4x4LookAt, podemos entender o uso da matriz transposta da seguinte maneira: Do lado esquerdo da equação teríamos uma transformação do vetor final, com relação à câmera, e do lado direito, uma ou mais transformações em relação ao sistema de coordenadas view. Como estamos interessados apenas no valor final de \vec{v'}, temos:

\displaystyle M_{lookup}\cdot\vec{v'}=M_{view}\cdot\vec{v}

\displaystyle \vec{v'} = M_{view}\cdot M^{T}_{lookup}\cdot\vec{v}

Pronto! Basta agora “afastar” a câmera de acordo com as coordenadas do vetor eye e a transformação está completa:

\displaystyle \vec{v'} = M_{translate}\cdot M_{view}\cdot M^{T}_{lookup}\cdot\vec{v}

Rotinas que foram retiradas do OpenGL 3+

O modo core profile do OpenGL retirou uma série de rotinas do seu repertório. Dentre elas, as rotinas de manipulação de matrizes (que foram colocadas à cargo dos shaders) e as rotinas da GLU (GL Utilities).

Lidar com matrizes é algo muito útil em computação gráfica. Permite uma maneira unificada e centralizada de realizar transformações geométricas, por exemplo. Como OpenGL “jogou fora” essas rotinas, provavelmente porque é mais fácil que o desenvolvedor use alguma biblioteca externa para evitar as limitações das pilhas de matrizes, temos que implementar nossas próprias. A primeira é a multiplicação de matrizes.

Vale lembrar que OpenGL lida, por padrão, com matrizes orientadas por coluna. Ou seja, diferente do padrão de C, onde matrizes são orientadas por linha (o incremento de um ponteiro que aponta para uma matriz apontará para o próximo item da linha). Uma maneira de visualizar a multiplicação de matrizes orientadas por coluna (column major) é este:

Coluna de A e linha de B, multiplicadas, dá item em C.

Coluna de A e linha de B, multiplicadas, dá item em C.

E, também, só para lembrar, a ordem da multiplicação é da direita para esquerda. Isso é assim porque a multiplicação de uma matriz por um vetor é feita assim:

\displaystyle \begin{bmatrix}x'\\y'\\z'\\w'\end{bmatrix}=\begin{bmatrix}m_1&m_5&m_9&m_{13}\\m_2&m_6&m_{10}&m_{14}\\m_3&m_7&m_{11}&m_{15}\\m_4&m_8&m_{12}&m_{16}\end{bmatrix}\cdot\begin{bmatrix}x\\y\\z\\w\end{bmatrix}

E matrizes podem ser encaradas como sendo um vetor de vetores. Cada coluna (mas matrizes column major) é um vetor.

A rotina de multiplicação de matrizes “padrão” é algo mais ou menos assim:

void Matrix4x4Multiply(float *mOut, const float *m1, const float *m2)
{
  int i, c, r;

  for (c = 0; c < 4; c++)
    for (r = 0; r < 4; r++)
    {
      mOut[c*4+r] = 0.0f;
      for (i = 0; i < 4; i++)
        mOut[c*4+r] += m1[c*4+i] * m2[i*4+r];
    }
}

Isso ai funciona muito bem, mas eis uma rotina um pouco mais veloz, usando SSE. Se você debugar a função abaixo, ela faz exatamente o que a função acima faz, mas usa 2 loops, ao invés de 3 – a rotina acima realiza 16 inicializações, 64 multiplicações e 64 adições. A rotina abaixo faz 16 adições e 16 multiplicações (de blocos de 4 floats):

void Matrix4x4Multiply(float *mOut, const float *m1, const float *m2)
{
  int i, j;
  __m128 a_column, b_column, r_column;

  for (i = 0; i < 16; i += 4)
  {
    b_column = _mm_set1_ps(b[i]);
    a_column = _mm_load_ps(a);
    r_column = _mm_mul_ps(a_column, b_column);

    for (j = 1; j < 4; j++)
    {
      b_column = _mm_set1_ps(b[i+j]);
      a_column = _mm_load_ps(&amp;a[j*4]);
      r_column = _mm_add_ps(_mm_mul_ps(a_column, b_column), r_column);
    }

    _mm_store_ps(&amp;out[i], r_column);
  }
}

Eis a rotina que substitui a glLoadIdentity(). Trata-se de uma simples cópia de 4 blocos de 4 floats, usando SSE:

void Matrix4x4Identity(float *m)
{
  static __m128 t[4] = { { 1.0f, 0.0f, 0.0f, 0.0f },
                         { 0.0f, 1.0f, 0.0f, 0.0f },
                         { 0.0f, 0.0f, 1.0f, 0.0f },
                         { 0.0f, 0.0f, 0.0f, 1.0f } };
  __m128 *p = (__m128 *)m;
  int i;

  for (i = 0; i < 4; i++)
    p[i] = t[i];
}

Uma função que não tem equivalente direto nas rotinas do OpenGL, mas é útil na manipulação de operações com matrizes, é a transposição. Transposição é uma transformação de uma matriz onde linhas viram colunas (e vice versa). A rotina pode ser implemtnada com dois loops, mas, usando SSE, e um pequeno macete, fica mais rápido:

void Matrix4x4Transpose(float *mOut, const float *mIn)
{
  __m128 m[4], t[4];

  m[0] = _mm_load_ps(mIn);
  m[1] = _mm_load_ps(mIn + 4);
  m[2] = _mm_load_ps(mIn + 8);
  m[3] = _mm_load_ps(mIn + 12);

  t[0] = _mm_unpacklo_ps(m[0], m[1]);
  t[1] = _mm_unpacklo_ps(m[2], m[3]);
  t[2] = _mm_unpackhi_ps(m[0], m[1]);
  t[3] = _mm_unpackhi_ps(m[2], m[3]);

  _mm_store_ps(mOut,      _mm_movelh_ps(t[0], t[1]));
  _mm_store_ps(mOut + 4,  _mm_movehl_ps(t[1], t[0]));
  _mm_store_ps(mOut + 8,  _mm_movelh_ps(t[2], t[3]));
  _mm_store_ps(mOut + 12, _mm_movehl_ps(t[3], t[2]));
}

As rotinas glFrustum, gluPerspective e gluLookAt também sumiram do OpenGL (do GL Utilities). Como todas as três são, de fato, manipulação de matrizes (as duas primeiras manipulam a matriz de projeção, a última, a matriz ModelView), ei-las:

/* Frustum projection Matrix. */
void Matrix4x4Frustum(float *mOut, float l, float r, float b, float t, float n, float f)
{
  float rml = r - l;
  float tmb = t - b;
  float fmn = f - n;
  float n2 = 2.0f * n;

  mOut[1] = mOut[2] = mOut[3] =
  mOut[4] = mOut[6] = mOut[7] =
  mOut[12] = mOut[13] = mOut[15] = 0.0f;

  mOut[0] = n2 / rml;

  mOut[5] = n2 / tmb;

  mOut[8] = (r+l) / rml;
  mOut[9] = (t+b) / tmb;
  mOut[10] = -(f+n) / fmn;
  mOut[11] = -1.0f;

  mOut[14] = -(f * n2) / fmn;
}

/* Substitui gluPerspective (cria um frustum 'simétrico'). */
void Matrix4x4Perspective(float *mOut, float fov, float aspectratio, float n, float f)
{
  float xmax, ymax;

  ymax = n * tanf(DEG2RAD(fov));
  xmax = ymax * aspectratio;
  Matrix4x4Frustum(mOut, -xmax, xmax, -ymax, ymax, n, f);
}

/* Substitui gluLookAt */
void Matrix4x4LookAt(float *mOut, const float *mIn, const float *eye, const float *target, const float *up)
{
  float fw[3], side[3], _up[3];
  float matrix2[16], resultMatrix[16];

  Vector3Sub(fw, target, eye);
  Vector3SelfNormalize(fw);

  /* o vetor 'side' aponta para a esquerda da câmera. */
  Vector3Cross(side, fw, up);
  Vector3SelfNormalize(side);

  /* _up é, necessariamente, unitário! E é recalculado aqui
     para evitar o inconveniente de termos um vetor 'up' não
     perpendicular aos vetores 'forward' e 'side'. */
  Vector3Cross(_up, side, fw);

  /* Poderia inserir os 3 vetores via função Vector3Copy, sobre uma
     matriz de identidade, mas assim parece ser mais rápido. */
  matrix2[0] = side[0];
  matrix2[4] = side[1];
  matrix2[8] = side[2];

  matrix2[1] = _up[0];
  matrix2[5] = _up[1];
  matrix2[9] = _up[2];

  /* Note que se a câmera está orientada para (0,-1,0), então 'target' terá
     um componnte z menor que o de 'eye'. O vetor 'forward' será negativo,
     por isso temos que invertê-lo. */
  matrix2[2] = -fw[0];
  matrix2[6] = -fw[1];
  matrix2[10] = -fw[2];

  matrix2[3] = matrix2[7] = matrix2[11] = 
  matrix2[12] = matrix2[13] = matrix2[14] = 0.0f;
  matrix2[15] = 1.0f;

  Matrix4x4Multiply(resultMatrix, mIn, matrix2);

  /* Aproveita 'matrix2' para criar a matriz de translação. */
  Matrix4x4Identity(matrix2);
  matrix2[3] = -eye[0];
  matrix2[7] = -eye[1];
  matrix2[11] = -eye[2];
  Matrix4x4Multiply(mOut, resultMatrix, matrix2);
}

A função Matrix4x4Frustum simplesmente implementa a matriz padrão do frustum genérico do OpenGL:

\displaystyle \begin{bmatrix}\frac{2n}{r-l} & 0 & \frac{r+l}{r-l} & 0 \\ 0 & \frac{2n}{t-b} & \frac{t+b}{t-b} & 0 \\ 0 & 0 & -\frac{f+n}{f-n} & -\frac{2fn}{f-n} \\ 0 & 0 & -1 & 0 \end{bmatrix}

A função Matrix4x4Perspective implementa exatamente o que gluPerspective faz: Ela cria um frustum simétrico, baseado no campo de visão (field of view), em graus, e no “aspect ratio” do plano “near” (4:3, 16:9, etc). Essa função calcula os valores do retângulo visível do plano near [os pontos (left, bottom) e (right, top)], onde a câmera aponta exatamente para o centro do retângulo, e chama Matrix4x4Frustum.

A função Matrix4x4LookAt é um pouco mais complicada… Ela calcula a translação para um novo sistema de coordenadas baseado nos 3 eixos onde a câmera está. Primeiro, calculamos o vetor que vai da posição da câmera (coordenada eye) para o alvo (target) e o chamamos de fw (de “forward”). Esse vetor é normalizado porque precisamos de um vetor unitário para compor o novo sistema de coordenadas. Neste ponto temos os vetores fwup e precisamos do terceiro vetor,, perpendicular a esses dois, que chamaremos de side. Por isso os produtos vetoriais. De posse dos 3 vetores, montamos a matriz transposta da matriz abaixo:

\displaystyle M_{lookAt}=\begin{bmatrix}s_x & up_x & -fw_x & 0 \\s_y & up_y & -fw_y & 0 \\s_z & up_z & -fw_z & 0 \\ 0 & 0 & 0 & 1\end{bmatrix}

Essa matriz (ou, melhor, sua transposta) “reorienta” o sistema de coordenadas de visualização. Precisamos afastar a câmera para a posição desejada e, por isso, uma matriz de translação final é usada:

\displaystyle T=\begin{bmatrix}1&0&0&-eye_x\\ 0&1&0&-eye_y\\ 0&0&1&-eye_z\\ 0&0&0&1\end{bmatrix}\cdot T_{view}\cdot\begin{bmatrix}s_x&s_y&s_z&0\\up_x&up_y&up_z&0\\-fw_x&-fw_y&-fw_z&0\\0&0&0&1\end{bmatrix}

A matriz T_{view}, inicial, geralmente é uma identidade.

Notou também que recalculamos o vetor \vec{up}? Isso torna a precaução de usar a dica da projeção de \vec{up} sobre \vec{forward} (aqui) desnecessária.

Usando papel e tesoura – recortando triângulos…

Imagine que você tenha um objeto complexo, feito com milhares de triângulos, por exemplo, este:

Mesh complexo de triângulos - formando um terreno.

Mesh complexo de triângulos – formando um terreno.

Parte dos triângulos no mesh não são visíveis para a câmera, como você pode observar. O que se deve fazer, neste caso, é recortar parte do mesh e desenhar apenas aqueles triângulos que são, de fato, visíveis. Este é um dos motivos da existência das rotinas de frustum culling. O frustum é uma figura trapezoide imaginária que existe na frente da câmera. Uma de seuas funções é usar os 6 planos que o definem para recortar os objetos que estão total ou parcialmente no seu interior e “jogar fora” todo o resto.

Uma nota de rodapé: Podemos usar sub-frustums, definidos por “portais” no mapa (que consta numa BSP) para determinar a visibilidade de objetos nas diversas regiões do mapa, a partir da posição da câmera, e recortar aqueles que não são visíveis, descartando-os.

Mas, ainda temos um problema: O que fazer com aqueles triângulos que são apenas parcialmente visíveis? Temos que pegar a tesoura e recortá-los para que caibam na parte vísivel, jogando fora as arestas! Neste processo de recorte alguns, senão muitos, triângulos que formam um objeto serão transformados em quadriláteros, como mostrado na figura abaixo:

Triângulo recortado por plano. Área cinza é descartada.

Triângulo recortado por plano. Área cinza é descartada.

O recorte é feito contra um dos planos do frustum. O que está dentro do frustum pode ser visto, o que está fora, não.

Na figura acima, ao ser recortado pelo plano o triângulo ABC é “transformado” no quadrilátero ABB’A’ (node o sentido anti-horário!). A mesma coisa não acontece se apenas um vértice do triângulo estiver do lado de “dentro” do frustum:

Triângulo recortao por plano. Área cinza é descatada.

Triângulo recortao por plano. Área cinza é descatada.

Sempre que temos apenas um dos vértices dentro do frustum vamos obter um triângulo, no recorte. Como regra geral, se apenas um vértice encontra-se do lado de fora do frustum, acabamos com um quadrilátero. Mas, como estamos lidando apenas com triângulos, precisaremos transformar um polígono (no caso, um quadrilátero) em um conjunto de triângulos. Esse processo chama-se tesselation. No nosso caso, fazer isso é simples. Basta escolher um dos vértices originais dentro do frustum e um dos novos vértices, no lado oposto, que está sobre o plano de clipping:

"Tesselando" o quadrilátero em 2 triângulos.

“Tesselando” o quadrilátero em 2 triângulos.

Outra coisa que você pode notar é que, não importa como os triângulos estejam dispostos, precisaremos obter dois novos vértices em dois dos lados (edges), exceto em um único caso: Se um dos vértices estiver precisamente sobre o plano de clipping. Neste caso, a área resultante também não precisará passar pelo processo de tesselation, já que a área “branca” continuará sendo um triângulo. A maneira de resolver essa possibilidade é considerar que se a distância do vértice ao plano for menor ou igual a zero, então ele está “atrás” do plano, deixando um vértice à frente o outro atrás.

Claro, que se todos os 3 vértices estiverem na frente do plano, nenhum clipping é necessário!

Bem… cadê o código para realizar o clippingtesselation? Você não vai encontrá-lo no código em C do projeto…. Isso é tarefa para os shaders. Se processarmos cada triângulo da cena sequencialmente, o tempo gasto será excessivo e os shaders são processos que executam em paralelo. No caso, os geometry shaderstesselation shaders operam em conjuntos de vértices. Se sua placa de vídeo pode usar 1024 shaders, por exemplo, cerca de 1024 triângulos podem ser clippadostesselados (se é que esses verbos existem) ao mesmo tempo.

Depois de explicar todos os tipos de shaders e seus funcionamentos, explicarei como implementar clipping.

Um raio interceptando um plano

Uma rotinazinha útil para ser usada com BSPs é a de interceptação de um raio com um plano. Um “raio” é definido aqui como um ponto no sistema de coordenadas associado com um vetor que dá uma direção. A equação é essa aqui:

\displaystyle \vec{r}=\dot{q}+t\vec{v}

Onde t é um valor usado para escalonar o vetor unitário \vec{v}. E \dot{q} é a posição inicial do raio.

Nós vimos, anteriormente, que um plano também segue uma equação:

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

Como queremos encontrar um ponto \dot{r} que está exatamente sobre a superfície do plano, substituímos a equação do raio na do plano e temos:

\displaystyle \begin{matrix}\vec{n}\cdot\dot{r}+d=0,\\\vec{n}\cdot(\dot{q}+t\vec{v})+d=0,\\\vec{n}\cdot\dot{q}+t(\vec{n}\cdot\vec{v})+d=0,\\t(\vec{n}\cdot\vec{v})=-(\vec{n}\cdot\dot{q}+d),\\t=\frac{-(\vec{n}\cdot\dot{q}+d)}{\vec{n}\cdot\vec{v}}\end{matrix}

Note que o numerador é a equação do plano, que equivale à chamada a função PlageGetDistanceFromPoint(). Assim, a equação para encontrar um ponto do raio de intercepta um plano é:

\displaystyle p_{intercept} = \dot{p}+\vec{v}\cdot\frac{-(\vec{n}\cdot\dot{q}+d)}{\vec{n}\cdot\vec{v}}

Onde \dot{p} é a coordenada inicial do raio, \vec{v} é o vetor unitário que dá a direção do raio e \vec{n} e ‘d’  fazem parte da equação do plano.

Isso nos dá o ponto p_{intercept}. Mas, existem dois detalhes: O produto escalar \vec{n}\cdot\vec{v}, se for zero, nos diz que o vetor direcional do raio é paralelo ao plano (chama-se vetor coplanar) e, portanto, não pode haver intersecção – este é o primeiro teste. O segundo teste é o numerador … se a distância do ponto \dot{q} até o plano for zero, significa que ele já está localizado no plano!

O primeiro teste faz todo sentido: Se um vetor qualquer for perpendicular a um vetor unitário, isto significa que o cosseno do ângulo entre eles é zero! Isso nos dá uma rotina útil para o módulo plane.c:

int PlaneInterceptedRay(float *vout, const float *plane, 
                        const float *pos, const float *vdir)
{
  float d, du;
  float u[3];

  /* Precisamos do vetor unitário de direção para calcular a escala. */
  Vector3Normalize(u, vdir);

  if ((du = Vector3Dot(plane, u)) == 0.0f)
    return -1;  /* coplanar? */

  /* Pertence ao plano? */
  if ((d = -PlaneGetDistanceFromPoint(plane, pos)) == 0.0f)
  {
    Vector3Copy(vout, pos);
    return 0;
  }

  Vector3Scale(vout, u, (d / du));
  Vector3Add(vout, pos);

  return 0;
}

Você pode estar se perguntando qual é a utilidade de descobrirmos um ponto de uma linha que intercepta um plano… Posso citar pelo menos dois: Recorte e visibilidade. Recorte é uma maneira de nos livrarmos de pedaços de triângulos que não são visíveis e quanto à visibilidade, diversos algorítmos podem usar isso: O desenho de sombras, por exemplo, depende da interceptação de feixes de luz com o plano onde a sombra aparecerá!

Á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.

Rotação em eixo arbitrário, SSE e Quaternions

Sem entrar em detalhes sobre a matemática dos quaternions (comece a ver aqui), eis a função de rotação de um vetor em torno de um eixo arbitrário (outro vetor), facilitada. Existe uma diferença com os códigos anteriores… agora uso SSE para lidar com os vetores. Os componentes do vetor são 4 floats: { x, y, z, w }:

__m128 vec4_rotateAroundAxis(__m128 v, __m128 axis, float angle)
{
  /* u = normalize(axis);
     vout = v*cos(angle) + (u x v)*sin(angle) + u*(u . v)*(1 - cos(angle)); */

  __m128 u;
  float c;

  c = cos(angle);
  u = vec4_normalize(axis);

  return _mm_add_ps(
    vec4_scale(v, c),
    _mm_add_ps(
      vec4_scale(vec4_cross(u, v), sin(angle)),
      vec4_scale(u, vec4_dot(u,v) * (1.0f * c))
    )
  );
}

As funções vec4_scale, vec4_cross, vec4_dot e vec4_normalize todas foram adaptadas das equivalentes vec3_*… Mas todas fazendo extensivo uso de SSE 4.1+. Em outro artigo coloco as funções aqui… Mas, você pode baixar o código com GIT:

$ git clone https://github.com/fredericopissarra/virtualworlds.git

Mais sobre o frustum…

Já falei sobre ele aqui. O frustum é um sólido trapezóide, uma pirâmide sem a parte de cima. O ponto focal (o topo da pirâmide) é onde está a câmera. Na projeção perspectiva essa pirâmide é calculada de acordo com o campo de visão (field of view), que é a amplitude angular que a câmera pode “enxergar”. Esse ângulo, no OpenGL, é fornecido como sendo o ângulo dos dois planos, superior e inferior, do frustum.

Ao chamar uma das funções, glFrustum ou gluPerspective, fornecemos os parâmetros para que o OpenGL calcule os 6 planos, baseados nas coordenadas do plano near e fare o campo de visão. O cálculo é simples e pode ser feito por semelhança de triangulos, dê uma olhada:

O frustum, com suas coordenadas calculadas por semelhança de triângulos.

Dá para entender agora porque as coordenadas (left,top,rright,bottom) fornecidas nas funçções glFrustum e gluPerspective são sempre relativas ao plano near . Com base neassas coordenadas e no field of view podemos calcular o plano far facilmente. Aqui vocẽ acha uma dedução do que essas funções fazem.

Com base nessa matriz podemos extrair o frustum de acordo com a orientação da câmera, levando em conta a matriz modelview.

/* plane.h */
...
typedef struct plane_s {
  float a, b, c, d;
} plane_t;

/* Veja observação após essa listagem */
void plane_normalize(plane_t *plane_p);

/* mat4.h */
...
#define MAT4_INDEX(col, row) ((row)*4 + (col))

void mat4_mult(float *mout_p, float *m1_p, float *m2_p);

/* frustum.c */

#include <GL/gl.h>
#include "mat4.h"
#include "plane.h"

void frustum_get_planes_from_opengl(plane_t *planes_p)
{
  float mat4_modelview[16];
  float mat4_projection[16];
  float mat4_temp[16];

  glGetFloatv( GL_MODELVIEW_MATRIX, mat4_modelview);
  glGetFloatv( GL_PROJECTION_MATRIX, mat4_projection);

  mat4_mult(mat4_temp, mat4_modelview, mat4_projection);

  /* left plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] + mat4_temp[MAT4_INDEX(0,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] + mat4_temp[MAT4_INDEX(0,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] + mat4_temp[MAT4_INDEX(0,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] + mat4_temp[MAT4_INDEX(0,3)];
  plane_normalize(planes_p++);

  /* right plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] - mat4_temp[MAT4_INDEX(0,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] - mat4_temp[MAT4_INDEX(0,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] - mat4_temp[MAT4_INDEX(0,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] - mat4_temp[MAT4_INDEX(0,3)];
  plane_normalize(planes_p++);

  /* bottom plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] + mat4_temp[MAT4_INDEX(1,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] + mat4_temp[MAT4_INDEX(1,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] + mat4_temp[MAT4_INDEX(1,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] + mat4_temp[MAT4_INDEX(1,3)];
  plane_normalize(planes_p++);

  /* top plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] - mat4_temp[MAT4_INDEX(1,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] - mat4_temp[MAT4_INDEX(1,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] - mat4_temp[MAT4_INDEX(1,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] - mat4_temp[MAT4_INDEX(1,3)];
  plane_normalize(planes_p++);

  /* near plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] + mat4_temp[MAT4_INDEX(2,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] + mat4_temp[MAT4_INDEX(2,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] + mat4_temp[MAT4_INDEX(2,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] + mat4_temp[MAT4_INDEX(2,3)];
  plane_normalize(planes_p++);

  /* far plane */
  planes_p->a = mat4_temp[MAT4_INDEX(3,0)] - mat4_temp[MAT4_INDEX(2,0)];
  planes_p->b = mat4_temp[MAT4_INDEX(3,1)] - mat4_temp[MAT4_INDEX(2,1)];
  planes_p->c = mat4_temp[MAT4_INDEX(3,2)] - mat4_temp[MAT4_INDEX(2,2)];
  planes_p->d = mat4_temp[MAT4_INDEX(3,3)] - mat4_temp[MAT4_INDEX(2,3)];
  plane_normalize(planes_p);
}

A função frustum_get_from_opengl() extrai os seis planos do frustum da combinação das matrizes, só que essas equações provavelmente não são normalizadas. Lembre-se que a equação do plano é:

Ax + By + Cz + D = 0

Onde (A,B,C) é um vetor unitário, normal (perpendicular) ao plano. Ao extrair os planos com a rotina acima esse vetor não será, provavelmente, unitário. Daí a necessidade de “normalizar” a equação:

void plane_normalize(plane_t *plane_p)
{
  float len;
  vec3_t *vnormal_p;

  vnormal_p = (vec3_t *)plane_p;
  len = vec3_length(vnormal_p);

  if (len != 0.0f)
  {
    plane_p->a /= len;
    plane_p->b /= len;
    plane_p->c /= len;
    plane_p->d /= len;
  }
}

Como a constante D é calculada a partir do vetor unitário, é necessário dividi-la pelo tamanho do vetor (A,B,C) original também.

Obter os planos do frustum é importante para fazermos testes de visibilidade. Um desses testes usa uma esfera onde o objeto está inscrito:

bool frustum_sphere_test(plane_t *planes_p, vec3_t *pt_p, float radius)
{
  int i;
  float distance;

  for (i = 0; i < 6; i++, planes_p++)
  {
    distance = plane_distance_from_point(planes_p, pt_p);

    /* Se o ponto está distante de -radius unidades
       do lado "negativo" do plano, então a esfera onde
       o objeto está inscrito está totalmente fora do frustum! */
    if (distance < -radius)
      return false;
  }

  /* A esfera está dentro ou parcialmente dentro do frustum. */
  return true;
}

O teste esférico é o primeiro teste de visibilidade que faremos (e o mais rápido).

Neste ponto você deve ter notado que o algoritmo tem um pequeno problema. Imagine uma esfera num dos vértices do frustum, cujo centro está distante do vértice exatamente do raio da esfera. Neste contexto a distância do centro do objeto (e da esfera onde está inscrito) está distante de dois planos por menos que o raio da esfera e o algoritmo dirá que o objeto está parcialmente dentro do frustum!

A solução para isso é considerar um delta_radius, uma fração do raio, somada ao raio da esfera — ou, visto de outra forma, como se o furstum tenha sido “estendido” por uma fração do raio da esfera. Essa solução eliminará o problema acima, mas criará outro. Objetos distantes de um dos planos pelo raio da esfera que estão, de fato, fora do frustum passarão a ser interpretados como estando parcialmente dentro dele. Assim, acabaríamos com mais objetos na lista de objetos possivelmente visíveis do que poderíamos querer.

Lembre-se que o teste esférico é apenas o primeiro teste. Usamos este teste para selecionar os objetos que potencialmente estão dentro do frustum para filtrá-los com testes mais precisos. No VirtualWorlds usarei testes com Object Oriented Bounding Boxes e Elipsóides, nesta ordem… Os testes com OOBB’s são mais precisos do que os esféricos, filtrando ainda mais os objetos possivelmente visíveis. E os testes com elipsoides filtrarão ainda mais. A vantagem é que esses outros testes terão que lidar com um subconjunto mais restrito de objetos…

Para facilitar, poderíamos mudar a rotina acima para marcar os objetos que são selecionados como potencialmente visíveis como estando totalmente ou parcialmente visíveis. Somente esses últimos seriam passados para a próxima fase de testes e assim por diante.