Photon mapping em uma cena simples. A técnica de photon mapping foi desenvolvida por Henrik W. Jensen(http://graphics.ucsd.edu/~henrik/). Este algoritmo é no código de Grant Schindler 2007.

http://www.cc.gatech.edu/~phlosoft/

Para compilar use: gcc -o cornel cornel.c -lglut

• Download cornell.c
• Enviar nova versão

Esconder código-fonte

/*
*   Copyright (C) 2008, 2008 Rafael Siqueira Telles Vieira
*
*   This program is free software; you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation; either version 2 of the License, or
*   (at your option) any later version.
*
*   This program is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*   See the GNU General Public License for more details.
*
*   License: http://www.lia.ufc.br/~rafaelstv/licenca_GPL_pt.txt
*
*   You should have received a copy of the GNU General Public License
*   along with this program; if not, write to the Free Software Foundation,
*   Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
*
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <GL/glut.h>

/** Cenario **/

int tamImagem = 400; // lagura e alturas
int nrTipos = 2; // tipos de objetos
int nrObjetos[] = {2,5}; // 2 esferas, 5 planos
float gOrigem[] = {0.0f,0.0f,0.0f};
float Luz[] = {0.0,1.2,3.75};
float esferas[][4] = {{1.0,0.0,4.0,0.5},{-0.6,1.0,4.5,0.5}}; // centro e raio
float planos[][2] = {{0, 1.5},{1, -1.5},{0, -1.5},{1, 1.5},{2,5.0}}; // eixo e distancia para origem

/** Globais do Photon Mapping **/

int nrPhotons = 1000;
int nrQuiques = 3;
float sqRadius = 0.70;
float exposicao = 50.0;
int numPhotons[][5] = {{0,0},{0,0,0,0,0}};
float photons[2][5][5000][3][3];

/** Globais do Raytracing **/

int gIntercepta = 0;
int gTipo;
int gId;
float gQdDist, gDist = -1.0;
float gPonto[] = {0.0, 0.0, 0.0};

/** Operações com Vetores **/

#define PROD_ESCALAR(vetor1, vetor2)\
   ((vetor1[0]*vetor2[0]) + (vetor1[1]*vetor2[1]) + (vetor1[2]*vetor2[2]))

#define MULT_ESCALAR(vetorS,vetorE, escalar)\
   vetorS[0] = vetorE[0] * escalar;\
   vetorS[1] = vetorE[1] * escalar;\
   vetorS[2] = vetorE[2] * escalar;

#define MULTIPLICA(resultado, vetor1, vetor2)\
   resultado[0] = vetor1[0] * vetor2[0];\
   resultado[1] = vetor1[1] * vetor2[1];\
   resultado[2] = vetor1[2] * vetor2[2];
 
#define SUBTRAI(resultado, vetor1, vetor2)\
   resultado[0] = vetor1[0] - vetor2[0];\
   resultado[1] = vetor1[1] - vetor2[1];\
   resultado[2] = vetor1[2] - vetor2[2];

#define SOMA(resultado, vetor1, vetor2)\
   resultado[0] = vetor1[0] + vetor2[0];\
   resultado[1] = vetor1[1] + vetor2[1];\
   resultado[2] = vetor1[2] + vetor2[2];

void RAND3(float vetor[])
{
   double random;
   random = -1.0 + ((double)rand()/(double)RAND_MAX)*2.0f;
   vetor[0] = random;
   random = -1.0 + ((double)rand()/(double)RAND_MAX)*2.0f;
   vetor[1] = random;
   random = -1.0 + ((double)rand()/(double)RAND_MAX)*2.0f;
   vetor[2] = random;
}

void NORMALIZA(float vetor[])
{
   float valor = sqrt(PROD_ESCALAR(vetor,vetor));
   vetor[0]/= valor;
   vetor[1]/= valor;
   vetor[2]/= valor;
}

int dentroRAIO(float a[], float b[], float sqraio)
{
   float c;
   float d = 0.0f;

   c = a[0] - b[0];
   d += c*c;
   if (d > sqraio) return 0;

   c = a[1] - b[1];
   d += c*c;
   if (d > sqraio) return 0;

   c = a[2] - b[2];
   d += c*c;
   if (d > sqraio) return 0;
   
   gQdDist = d;
   return 1;
}

/** Intersecoes **/

void verificaDistancia(float Dist, int tipo, int id)
{
   if (Dist < gDist && Dist > 0.0)
   {
      gTipo = tipo;
      gDist = Dist;
      gId = id;
      gIntercepta = 1;
   }
}

void raioEsfera(int id, float direcao[], float origem[])
{
   float t[3]={0.0,0.0,0.0}; 
   float raio = esferas[id][3];

   SUBTRAI(t,esferas[id],origem);

   float A = PROD_ESCALAR(direcao,direcao);
   float B = -2.0f*PROD_ESCALAR(direcao,t);
   float C = PROD_ESCALAR(t,t) - (raio*raio);
   float D = B*B - 4*A*C;
   
   if (D > 0.0)
   {
      float sinal = (C < -0.0001) ? 1: -1;
      float Dist1 = (-B + sinal*sqrt(D))/(2*A);   
      verificaDistancia(Dist1,0,id);
   }
}

void raioPlano(int id, float direcao[], float origem[])
{
   int eixo = (int) planos[id][0];
   if (direcao[eixo] != 0.0)
   {
      float Dist = (planos[id][1] - origem[eixo]) / direcao[eixo];
      verificaDistancia(Dist,1,id);
   } 
}

void raioObjeto(int tipo, int id, float direcao[], float origem[])
{
   if (tipo == 0) 
      raioEsfera(id,direcao,origem);
   else
      raioPlano(id,direcao,origem);
}


/** Normal **/

void esferaNormal(float normal[], int id, float P[])
{
   SUBTRAI(normal,P,esferas[id]);
   NORMALIZA(normal);
}

void planoNormal(float normal[], int id, float P[], float origem[])
{
   int eixo = (int) planos[id][0];
   normal[0] = 0.0;
   normal[1] = 0.0;
   normal[2] = 0.0;
   normal[eixo] = origem[eixo] - planos[id][1];
   NORMALIZA(normal);
}

void superficieNormal(float normal[], int tipo, int id, float P[], float dentro[])
{
   if (tipo == 0)
      esferaNormal(normal,id,P);
   else
      planoNormal(normal,id,P,dentro);
}

/** Raytracing **/

void tracarRaio(float raio[], float origem [])
{
   int t,i;
   gIntercepta = 0;
   gDist = 9999.9;
   for (t=0; t<nrTipos;t++)
      for (i=0; i<nrObjetos[t];i++)
         raioObjeto(t,i,raio,origem);

}

void reflita(float rRefletido[], float rIncidente[], float pontoPartida[])
{
   float normal[3]={0.0,0.0,0.0};
   float resultado[3]={0.0,0.0,0.0};
   superficieNormal(normal,gTipo,gId,gPonto,pontoPartida);
   MULT_ESCALAR(resultado,normal,PROD_ESCALAR(normal,rIncidente)*2.0);
   SUBTRAI(rRefletido,rIncidente,resultado);
   NORMALIZA(rRefletido);
}

/** Photon Mapping **/

void armazenaPhoton(int tipo, int id, float localizacao[], float direcao[], float energia[])
{
   photons[tipo][id][numPhotons[tipo][id]][0][0] = localizacao[0];
   photons[tipo][id][numPhotons[tipo][id]][0][1] = localizacao[1];
   photons[tipo][id][numPhotons[tipo][id]][0][2] = localizacao[2];
   photons[tipo][id][numPhotons[tipo][id]][1][0] = direcao[0];
   photons[tipo][id][numPhotons[tipo][id]][1][1] = direcao[1];
   photons[tipo][id][numPhotons[tipo][id]][1][2] = direcao[2];
   photons[tipo][id][numPhotons[tipo][id]][2][0] = energia[0];
   photons[tipo][id][numPhotons[tipo][id]][2][1] = energia[1];
   photons[tipo][id][numPhotons[tipo][id]][2][2] = energia[2];
   numPhotons[tipo][id]++;
}

void filtreCor(float corS[], float corE[], float r, float g, float b)
{   
   int c;
   corS[0] = r; 
   corS[1] = g;
   corS[2] = b;
   for (c = 0; c < 3; c++)
      if (corS[c] > corE[c])
      corS[c] = corE[c]; 
}
void retCor(float corS[], float corE[], int tipo, int id)
{
   if (tipo == 1 && id == 0)    { filtreCor(corS, corE, 0.0, 1.0, 0.0); }
   else if (tipo == 1 && id == 2)  { filtreCor(corS, corE, 1.0, 0.0, 0.0); }
   else            { filtreCor(corS, corE, 1.0, 1.0, 1.0); }
}

void reunaPhotons(float cor[], float p[], int tipo, int id)
{
   int i;
   cor[0] = 0.0;
   cor[1] = 0.0;
   cor[2] = 0.0;
   float N[3];
   float resultado[3];

   superficieNormal(N,tipo,id,p,gOrigem);
   for (i = 0; i < numPhotons[tipo][id]; i++)
   {
      if (dentroRAIO(p,photons[tipo][id][i][0],sqRadius)) 
      {
         
         float max = PROD_ESCALAR(N,photons[tipo][id][i][1])*-1.0f;
         if (max < 0.0f)
            max = 0.0f;
         max *= ((1.0 - sqrt(gQdDist))/exposicao);
         resultado[0] = 0.0f;
         resultado[1] = 0.0f;
         resultado[2] = 0.0f;

         MULT_ESCALAR(resultado,photons[tipo][id][i][2],max);
         SOMA(cor, cor, resultado);          
      }
   }
}

void shadowPhoton(float raio[])
{
   float shadow[] = {-0.25,-0.25,-0.25};
   float ponto[] = {gPonto[0], gPonto[1], gPonto[2]};
   int tipo = gTipo, id = gId, dist = gDist;
   float pontoMovido[3]={0.0,0.0,0.0},resultado[3]={0.0,0.0,0.0};
   MULT_ESCALAR(resultado,raio,0.001);
   SOMA(pontoMovido,resultado,gPonto);
   tracarRaio(raio,pontoMovido);
   float pontoSombra[3]={0.0,0.0,0.0};
   resultado[0] = 0.0f;
   resultado[1] = 0.0f;
   resultado[2] = 0.0f;
   MULT_ESCALAR(resultado,raio,gDist);
    SOMA(pontoSombra, resultado,pontoMovido);
   armazenaPhoton(gTipo,gId,pontoSombra,raio,shadow);
   gPonto[0]=ponto[0];
   gPonto[1]=ponto[1];
   gPonto[2]=ponto[2];
   gTipo = tipo;
   gId = id;   
   gDist = dist;
}

/** colorir **/

void computeCorPixel(float cor[], float x, float y)
{
   float raio[3] = {(x/tamImagem)-0.5, -((y/tamImagem)-0.5),1.0};
   int tipo;
   int id;
   tracarRaio(raio,gOrigem);
   if (gIntercepta)
   {
      MULT_ESCALAR(gPonto,raio,gDist);
      if (gTipo == 0 && gId == 1)
      {
         float refletido[3]={0.0f,0.0f,0.0f};
         reflita(refletido,raio,gOrigem);
         tracarRaio(refletido,gPonto);
         if (gIntercepta)
         {
            float resultado[3]={0.0f,0.0f,0.0f};
            MULT_ESCALAR(resultado,refletido,gDist);
            SOMA(gPonto,resultado,gPonto);
         }
      }
      reunaPhotons(cor,gPonto,gTipo,gId);
   }
}

void emitaPhotons()
{
   int i,t;
   srand(time(0));
   for (t = 0; t<nrTipos; t++)
      for (i=0;i<nrObjetos[t];i++)
         numPhotons[t][i] = 0;

   for (i=0;i<nrPhotons;i++)
   {
      int quiques = 1;
      float cor[] = {1.0,1.0,1.0};
      float raio[3];

      RAND3(raio);
      NORMALIZA(raio);
      float pontoAnt[] = { Luz[0], Luz[1], Luz[2] };

      while (pontoAnt[1] >= Luz[1])
      {
         float resultado[3];
         RAND3(resultado);
         NORMALIZA(resultado);
         MULT_ESCALAR(resultado,resultado,0.75);
         SOMA(pontoAnt,Luz,resultado);
      }

      if (abs(pontoAnt[0]) > 1.5 || (abs(pontoAnt[1]) > 1.2) 
      || dentroRAIO(pontoAnt, esferas[0],esferas[0][3]*esferas[0][3]))
         quiques = nrQuiques + 1;

      tracarRaio(raio,pontoAnt);
         
      float resultado[3];
      while (gIntercepta && quiques <= nrQuiques)
      {
         resultado[0] = 0.0f;
         resultado[1] = 0.0f;
         resultado[2] = 0.0f;
         MULT_ESCALAR(resultado,raio,gDist);
         SOMA(gPonto,resultado,pontoAnt);
         retCor(cor,cor,gTipo,gId);
         MULT_ESCALAR(cor,cor,1.0/sqrt(quiques));
         armazenaPhoton(gTipo,gId,gPonto,raio,cor);
         shadowPhoton(raio);
         float raioR[3];
         reflita(raioR, raio, pontoAnt);
         tracarRaio(raioR,gPonto);
         pontoAnt[0] = gPonto[0];
         pontoAnt[1] = gPonto[1];
         pontoAnt[2] = gPonto[2];
         quiques++;
      }
   }
}

void Inicia()
{
   emitaPhotons();   
}

/** OpenGL **/
void desenha(void)
{
   glMatrixMode(GL_PROJECTION);
       glLoadIdentity();
       gluOrtho2D(0,tamImagem,0,tamImagem);

   glClearColor(1.0f,1.0f,1.0f,1.0f);
    glClear(GL_COLOR_BUFFER_BIT);
   int x,y;
   float cor[3];
   for (x=0;x<tamImagem;x++)
      for (y=0;y<tamImagem;y++)
      {
         cor[0]=0.0f;
         cor[1]=0.0f;
         cor[2]=0.0f;
         computeCorPixel(cor,x,y);
         glRasterPos2f(x,y); 
          glDrawPixels(1,1,GL_RGB,GL_FLOAT,cor);
      }
    glFlush;
}

int main(int argc, char** argv)
{
    glutInit(&argc,argv);
   glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
   glutInitWindowSize(tamImagem,tamImagem);
   glutInitWindowPosition(0,0);
   glutCreateWindow("Photon Mapping");
   glutDisplayFunc(desenha);
   Inicia();
   glutMainLoop();

   return 0;
}

Algoritmo de Ordenação Radix

Pilha

[C] Qsort - Ordenação

Potências de 2

merge sort