Thanks, my friend. I don’t know if my allocation is correct. I’m stuck with that problem since 5 A.M of thursday, rsrsrs.
If you guys could help me figure it out. Here is my code, I’m realy sorry for it be kind messy.
/* *
* Copyright 1993-2012 NVIDIA Corporation. All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "rply.h"
/**
* This macro checks return value of the CUDA runtime call and exits
* the application if the call failed.
*/
#define CUDA_CHECK_RETURN(value) { \
cudaError_t _m_cudaStat = value; \
if (_m_cudaStat != cudaSuccess) { \
fprintf(stderr, "Error %s at line %d in file %s\n", \
cudaGetErrorString(_m_cudaStat), __LINE__, __FILE__); \
exit(1); \
} }
__host__ __device__ unsigned int bitreverse(unsigned int number) {
number = ((0xf0f0f0f0 & number) >> 4) | ((0x0f0f0f0f & number) << 4);
number = ((0xcccccccc & number) >> 2) | ((0x33333333 & number) << 2);
number = ((0xaaaaaaaa & number) >> 1) | ((0x55555555 & number) << 1);
return number;
}
/**
* CUDA kernel function that reverses the order of bits in each element of the array.
*/
//MACROS
#define FLOAT double
#define QTD_LUZ 1
#define HEIGH 600
#define WIDTH 800
#define TAM_RAIO 100000
#define LIMITE 5
#define ERR_MARGIN 1e-6
#define MIN_TRIANGULOS 10
#define PIx180 0.017453293
typedef struct obj3d{
FLOAT x,y,z;
} obj3d;
typedef struct color{
FLOAT r,g,b;
} color;
typedef struct ipoint{
obj3d pos;
obj3d vrefl;
obj3d normal;
FLOAT dist;
long id;
int id_tri;
} ipoint;
typedef struct triang{
obj3d A;
obj3d B;
obj3d C;
obj3d normal;
obj3d baricentro;
color col;
FLOAT dif;
FLOAT spec;
FLOAT shi;
FLOAT a;
FLOAT refl;
} triang;
typedef struct esf{
obj3d centro;
FLOAT r;
color col;
FLOAT dif;
FLOAT spec;
FLOAT shi;
FLOAT a;
FLOAT refl;
triang lista_triang[MIN_TRIANGULOS];
long qtd_triangulos;
} esf;
typedef struct camera{
obj3d eye, dir, up;
FLOAT fov;
} camera;
typedef struct ray{
obj3d orig, dir;
} ray;
typedef struct light{
obj3d pos;
color dif;
color spec;
} light;
__host__ __device__ obj3d cross_product(obj3d v1, obj3d v2) {
obj3d res;
res.x = v1.y * v2.z - v1.z * v2.y;
res.y = v1.z * v2.x - v1.x * v2.z;
res.z = v1.x * v2.y - v1.y * v2.x;
return res;
}
//MACROS: Operações geométricas básicas
#define DET(A, B, C) ((((A).x * (B).y * (C).z) + ((B).x * (C).y * (A).z) + ((C).x * (A).y * (B).z)) - (((C).x * (B).y * (A).z) + ((A).x * (C).y * (B).z) + ((B).x * (A).y * (C).z)))
#define SQ(x) ((x) * (x))
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define DOT(a, b) ((a).x * (b).x + (a).y * (b).y + (a).z * (b).z)
#define NORMALIZE(a) do {\
FLOAT len = sqrt(DOT(a, a));\
(a).x /= len; (a).y /= len; (a).z /= len;\
} while(0);
__device__ color trace(ray raio, int lim, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes);
__device__ color shade_esf(esf esfera, ipoint S1, int lim, ray raio, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes);
__device__ color shade(triang tri, ipoint S1, int lim, ray raio, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes);
FLOAT aspect = 1.333333;
color *h_imagem, *imagem;
obj3d *vertices;
triang *faces;
esf *d_esferas, *h_esferas;
esf *d_esf_cab, *h_esf_cab;
long jind_v=0, jind=0, vert_aux1, vert_aux2, vert_aux3;
long qtd_nodes;
int altura;
__host__ __device__ obj3d rotacao(obj3d aux, int tipo, FLOAT t){
obj3d A;
t*=PIx180;
switch(tipo)
{ case 0:
A=aux;
break;
case 1:
(A).y = (((aux).y * (cos(t))) - ((aux).z * (sin(t))));
(A).z = (((aux).y * (sin(t))) + ((aux).z * (cos(t))));
break;
case 2:
(A).x = (((aux).x * (cos(t))) + ((aux).z * (sin(t))));
(A).z = (-((aux).x * (sin(t))) + ((aux).z * (cos(t))));
break;
case 3:
(A).x = (((aux).x * (cos(t))) - ((aux).y * (sin(t))));
(A).y = (((aux).x * (sin(t))) + ((aux).y * (cos(t))));
break;
}
return A;
}
__host__ __device__ obj3d escala(obj3d aux, FLOAT sx, FLOAT sy, FLOAT sz){
obj3d A;
(A).x = (aux).x * (sx);
(A).y = (aux).y * (sy);
(A).z = (aux).z * (sz);
return A;
}
__host__ __device__ obj3d trans(obj3d aux, FLOAT tx, FLOAT ty, FLOAT tz){
obj3d A;
(A).x = (aux).x + (tx);
(A).y = (aux).y + (ty);
(A).z = (aux).z + (tz);
return A;
}
triang construc_triang(obj3d A, obj3d B, obj3d C,
FLOAT r, FLOAT g, FLOAT b, FLOAT a,
FLOAT spec, FLOAT dif, FLOAT shi, FLOAT refl){
triang tri;
obj3d tri_BA;
obj3d tri_CA;
tri.A.x=A.x; tri.A.y=A.y; tri.A.z=A.z;
tri.B.x=B.x; tri.B.y=B.y; tri.B.z=B.z;
tri.C.x=C.x; tri.C.y=C.y; tri.C.z=C.z;
tri.baricentro.x= (A.x + B.x + C.x)/3;
tri.baricentro.y= (A.y + B.y + C.y)/3;
tri.baricentro.z= (A.z + B.z + C.z)/3;
tri.col.r=r; tri.col.g=g; tri.col.b=b; tri.a=a; tri.refl=refl;
tri.spec=spec;
tri.dif=dif;
tri.shi=shi;
tri_BA.x=tri.B.x - tri.A.x;
tri_BA.y=tri.B.y - tri.A.y;
tri_BA.z=tri.B.z - tri.A.z;
tri_CA.x=tri.C.x - tri.A.x;
tri_CA.y=tri.C.y - tri.A.y;
tri_CA.z=tri.C.z - tri.A.z;
tri.normal=cross_product(tri_BA, tri_CA);
NORMALIZE(tri.normal);
return tri;
}
esf construct_esf(FLOAT x, FLOAT y, FLOAT z, FLOAT raio, FLOAT r, FLOAT g, FLOAT b, FLOAT a, FLOAT spec, FLOAT dif, FLOAT shi, FLOAT refl){
esf esf_aux;
esf_aux.centro.x=x;
esf_aux.centro.y=y;
esf_aux.centro.z=z;
esf_aux.r=raio;
esf_aux.col.r=r;
esf_aux.col.g=g;
esf_aux.col.b=b;
esf_aux.a=a;
esf_aux.spec=spec;
esf_aux.dif=dif;
esf_aux.shi=shi;
esf_aux.refl=refl;
return esf_aux;
}
void dump_triang(triang tri, int i){
printf("triangulo %d:\n",i);
printf("Vertices A(%f,%f,%f) B(%f,%f,%f) C(%f,%f,%f)\n", tri.A.x, tri.A.y, tri.A.z,
tri.B.x, tri.B.y, tri.B.z,
tri.C.x, tri.C.y, tri.C.z);
printf("Normal N(%f,%f,%f)\n",tri.normal.x,tri.normal.y,tri.normal.z);
printf("Baricentro N(%f,%f,%f)\n",tri.baricentro.x,tri.baricentro.y,tri.baricentro.z);
/*printf("componentes rgba(%f,%f,%f,%f) spec=%f shi=%f dif=%f\n\n", tri.col.r, tri.col.g, tri.col.b, tri.a,
tri.spec, tri.shi, tri.dif);*/
}
void dump_camera(FLOAT m[12], camera cam){
int i;
printf("Camera:\n");
printf("Posicao (%f, %f, %f) Direcao(%f, %f, %f) Up(%f, %f, %f) fov=%f\n\n", cam.eye.x, cam.eye.y, cam.eye.z,
cam.dir.x, cam.dir.y, cam.dir.z,
cam.up.x, cam.up.y, cam.up.z,
cam.fov);
printf("Matriz Transformação\n");
for(i=0;i<9;i+=3){
printf("%f %f %f\n", m[i], m[i+1], m[i+2]);
}
}
void dump_ray(ray raio, int i){
printf("Raio %d\n", i);
printf("O(%f, %f, %f) Dir(%f, %f, %f)\n",raio.orig.x,raio.orig.y, raio.orig.z,
raio.dir.x, raio.dir.y, raio.dir.z);
}
/* calculate reflection vector */
__device__ obj3d get_vrefletido(obj3d v, obj3d n) {
obj3d refl;
FLOAT dot = v.x * n.x + v.y * n.y + v.z * n.z;
refl.x = -(2.0 * dot * n.x - v.x);
refl.y = -(2.0 * dot * n.y - v.y);
refl.z = -(2.0 * dot * n.z - v.z);
return refl;
}
/*
* Função criada para verificar se um raio intercepta um triângulo.
* Recebe como parâmetros o raio, o triangulo em questão e um ipoint
* como referência que armazenará os dados do ponto interceptado.
* A função retorna 1 se ocorreu interceptação e 0 do contrário.
*/
__device__ int intercepta(ray raio, triang tri, ipoint *S){
//intercepta plano?
obj3d aux1, aux2, P, u, v, w;
FLOAT dot1, dot2, dot3, dot4, dot5, dot6, r, s, t;
aux1.x=tri.A.x - raio.orig.x;//V0 - P0
aux1.y=tri.A.y - raio.orig.y;
aux1.z=tri.A.z - raio.orig.z;
aux2.x=raio.dir.x - raio.orig.x;//P1 - P0
aux2.y=raio.dir.y - raio.orig.y;
aux2.z=raio.dir.z - raio.orig.z;
dot2=DOT(tri.normal,aux2);
//Não intercepta plano, muito menos triângulo.
if(dot2==0)
return 0;
dot1=DOT(tri.normal,aux1);
r=dot1/dot2;
//Não intercepta plano, muito menos triângulo.
if(r<0)
return 0;
P.x = raio.orig.x + r * (aux2.x);
P.y = raio.orig.y + r * (aux2.y);
P.z = raio.orig.z + r * (aux2.z);
w.x=P.x - tri.A.x;
w.y=P.y - tri.A.y;
w.z=P.z - tri.A.z;
u.x=tri.B.x - tri.A.x;//v1-v0
u.y=tri.B.y - tri.A.y;
u.z=tri.B.z - tri.A.z;
v.x=tri.C.x - tri.A.x; //v2-v0
v.y=tri.C.y - tri.A.y;
v.z=tri.C.z - tri.A.z;
dot1= DOT(u,v);
dot2= DOT(w,v);
dot3= DOT(v,v);
dot4= DOT(w,u);
dot5= DOT(u,u);
dot6= (dot1*dot1) - (dot5 * dot3);
s = ((dot1 * dot2) - (dot3 * dot4)) / dot6;
t = ((dot1 * dot4) - (dot5 * dot2)) / dot6;
//Reaproveitando variáveis
u.x=tri.B.x - P.x;
u.y=tri.B.y - P.y;
u.z=tri.B.z - P.z;
v.x=tri.C.x - P.x;
v.y=tri.C.y - P.y;
v.z=tri.C.z - P.z;
w.x=P.x - raio.orig.x;
w.y=P.y - raio.orig.y;
w.z=P.z - raio.orig.z;
S->dist=sqrt(DOT(w, w));
S->normal=cross_product(u, v);
if(DOT(raio.dir,S->normal)>0){
S->normal.x*=-1;
S->normal.y*=-1;
S->normal.z*=-1;
}
//Intercepta Triângulo
if(s>=0 && t>=0 && ((s + t) <= 1) && S->dist>ERR_MARGIN){
NORMALIZE(S->normal);
NORMALIZE(raio.dir);
S->vrefl=get_vrefletido(raio.dir,S->normal);
NORMALIZE(S->vrefl);
S->pos=P;
return 1;
}
return 0;
}
__device__ int intercepta_esf(ray raio, esf esfera, ipoint *sp){
FLOAT a, b, c, d, sqrt_d, t1, t2;
a = SQ(raio.dir.x) + SQ(raio.dir.y) + SQ(raio.dir.z);
b = 2.0 * raio.dir.x * (raio.orig.x - esfera.centro.x) +
2.0 * raio.dir.y * (raio.orig.y - esfera.centro.y) +
2.0 * raio.dir.z * (raio.orig.z - esfera.centro.z);
c = SQ(esfera.centro.x) + SQ(esfera.centro.y) + SQ(esfera.centro.z) +
SQ(raio.orig.x) + SQ(raio.orig.y) + SQ(raio.orig.z) +
2.0 * (-esfera.centro.x * raio.orig.x - esfera.centro.y * raio.orig.y - esfera.centro.z * raio.orig.z) - SQ(esfera.r);
if((d = SQ(b) - 4.0 * a * c) < 0.0) return 0;
sqrt_d = sqrt(d);
t1 = (-b + sqrt_d) / (2.0 * a);
t2 = (-b - sqrt_d) / (2.0 * a);
if((t1 < ERR_MARGIN && t2 < ERR_MARGIN) || (t1 > 1.0 && t2 > 1.0)) return 0;
if(sp) {
if(t1 < ERR_MARGIN) t1 = t2;
if(t2 < ERR_MARGIN) t2 = t1;
sp->dist = t1 < t2 ? t1 : t2;
sp->pos.x = raio.orig.x + raio.dir.x * sp->dist;
sp->pos.y = raio.orig.y + raio.dir.y * sp->dist;
sp->pos.z = raio.orig.z + raio.dir.z * sp->dist;
sp->normal.x = (sp->pos.x - esfera.centro.x) / esfera.r;
sp->normal.y = (sp->pos.y - esfera.centro.y) / esfera.r;
sp->normal.z = (sp->pos.z - esfera.centro.z) / esfera.r;
sp->vrefl = get_vrefletido(raio.dir, sp->normal);
NORMALIZE(sp->vrefl);
}
return 1;
}
/*
* Função utilizada para calcular a posição do raio no sistema de coordenadas
* do view frame. Serve tambpem para calcular sua direção. Recebe como parâmetro
* o valor da linha e coluna corrente da imagem a ser criada, e a quantidade de raios
* por pixel Por fim, retorna uma estrutura do tipo raio com sua posição e direção
* calculados, o qual será utilizado na função trace.
*/
__device__ ray get_p_ray(int x, int y, int ray_per_pixel, FLOAT m[12]){
ray raio;
obj3d dir_aux;
raio.dir.x = ((FLOAT) WIDTH) * (((FLOAT) x / (FLOAT) WIDTH) - 0.5);
raio.dir.y = -(((FLOAT) HEIGH)* (((FLOAT) y / (FLOAT) HEIGH) - 0.5)) ;
//raio.dir.z = -1;
dir_aux.x = raio.dir.x * m[0] + raio.dir.y * m[3] - m[6];
dir_aux.y = raio.dir.x * m[1] + raio.dir.y * m[4] - m[7];
dir_aux.z = raio.dir.x * m[2] + raio.dir.y * m[5] - m[8];
dir_aux.x-=m[9];
dir_aux.y-=m[10];
dir_aux.z-=m[11];
NORMALIZE(dir_aux);
dir_aux.x*=TAM_RAIO;
dir_aux.y*=TAM_RAIO;
dir_aux.z*=TAM_RAIO;
raio.orig.x=m[9];
raio.orig.y=m[10];
raio.orig.z=m[11];
raio.dir.x=dir_aux.x + raio.orig.x;
raio.dir.y=dir_aux.y + raio.orig.y;
raio.dir.z=dir_aux.z + raio.orig.z;
return raio;
}
/*
* Função para posicionar a câmera. Está receberá os dados necessários
* para posicionar a câmera: eye, traget, up e fov. Com os dados
* eye, target e up calculará o view coordinate frame. Por fim, cria
* uma matriz, a qual é uma variável global, que servirá para posteriormente
* expressar a posição inicial do raio em termos desse sistema de coordenadas.
*/
__device__ void posiciona_camera(FLOAT eye_x, FLOAT eye_y, FLOAT eye_z,
FLOAT targ_x, FLOAT targ_y, FLOAT targ_z,
FLOAT up_x, FLOAT up_y, FLOAT up_z,
FLOAT fov, FLOAT m[12]){
camera cam;
obj3d Ortho_cam_i, Ortho_cam_k, Ortho_cam_j;
cam.eye.x=eye_x; cam.eye.y=eye_y; cam.eye.z=eye_z;
cam.dir.x=targ_x - eye_x;
cam.dir.y=targ_y - eye_y;
cam.dir.z=targ_z - eye_z;
cam.up.x=up_x; cam.up.y=up_y; cam.up.z=up_z;
cam.fov=fov;
Ortho_cam_k = cam.dir;
NORMALIZE(Ortho_cam_k);
Ortho_cam_k.x*=-1; Ortho_cam_k.y*=-1; Ortho_cam_k.z*=-1;
Ortho_cam_i = cross_product(cam.up, Ortho_cam_k);
Ortho_cam_j = cross_product(Ortho_cam_k, Ortho_cam_i);
m[0] = Ortho_cam_i.x; m[3] = Ortho_cam_j.x; m[6] = Ortho_cam_k.x;
m[1] = Ortho_cam_i.y; m[4] = Ortho_cam_j.y; m[7] = Ortho_cam_k.y;
m[2] = Ortho_cam_i.z; m[5] = Ortho_cam_j.z; m[8] = Ortho_cam_k.z;
m[9] = cam.eye.x; m[10] = cam.eye.y; m[11] = cam.eye.z;
}
__device__ void cria_luz(int i, FLOAT dif_r, FLOAT dif_g, FLOAT dif_b,
FLOAT spec_r, FLOAT spec_g, FLOAT spec_b,
FLOAT x, FLOAT y, FLOAT z, light luz[QTD_LUZ]){
luz[i].dif.r=dif_r; luz[i].dif.g=dif_g; luz[i].dif.b=dif_b;
luz[i].spec.r=spec_r; luz[i].spec.g=spec_g; luz[i].spec.b=spec_b;
luz[i].pos.x=x; luz[i].pos.y=y; luz[i].pos.z=z;
}
void create_image_ppm(char *f_name){
int i, j;
FILE *fp = fopen(f_name, "wb"); /* b - binary mode */
(void) fprintf(fp, "P6\n%d %d\n255\n", WIDTH, HEIGH);
static unsigned char cor[3];
for (i = 0; i < HEIGH; i++)
{
for (j = 0; j < WIDTH; j++)
{
cor[0] = h_imagem[i * WIDTH + j].r; /* red */
cor[1] = h_imagem[i * WIDTH + j].g; /* green */
cor[2] = h_imagem[i * WIDTH + j].b; /* blue */
(void) fwrite(cor, 1, 3, fp);
}
}
(void) fclose(fp);
}
/*
* Função criada para caminhar recursivamente a spheere-tree até encontrar
* o nó que possua a quantidade de triângulos definida como mínima. Depois, procura
* sequencialmente cada triângulo, e armazena os dados referentes ao que se encontrar
* mais próximo à origem do raio.
* Recebe como parâmetro o id da raiz e o raio. Retorna uma estrutura do tipo ipoint que
* contém dados sobre o ponto do objeto onde ocorreu a interseção.
*/
__device__ ipoint get_esf_node(ray raio, long id_no, esf *esf_cab, long qtd_nodes){
long i_dir=(id_no+1) * 2; long k;
long i_esq=(id_no+1) * 2 - 1;
ipoint S1; S1.dist=TAM_RAIO+100; S1.id=-1;
ipoint S2; S2.dist=TAM_RAIO+100; S2.id=-1;
ipoint dummy;
if(id_no>qtd_nodes){
return S2;
}
if(intercepta_esf(raio, esf_cab[id_no],&dummy)==1){
if(esf_cab[id_no].qtd_triangulos <= MIN_TRIANGULOS){
//printf("entreoi id_no=%ld\n",id_no);
for(k=0;k<esf_cab[id_no].qtd_triangulos; k++){ //Procura triângulos.
if(intercepta(raio, esf_cab[id_no].lista_triang[k] ,&S1)==1 && S1.dist<S2.dist){
S2=S1;
S2.id=id_no;
S2.id_tri=k;
}
}
return S2;
}
else{
S1=get_esf_node(raio,i_esq,esf_cab, qtd_nodes);
S2=get_esf_node(raio,i_dir,esf_cab, qtd_nodes);
}
}
else
return S2;
if(S1.dist <= S2.dist)
return S1;
else
return S2;
}
/*__device__ ipoint get_esf_node_simpl(ray raio, long id_no, esf *esf_cab){
long i_dir=(id_no+1) * 2; long k;
long i_esq=(id_no+1) * 2 - 1;
ipoint S1; S1.dist=TAM_RAIO+100; S1.id=-1;
ipoint S2; S2.dist=TAM_RAIO+100; S2.id=-1;
ipoint dummy;
/*if(i_dir>qtd_nodes)
return S2;
if(i_esq>qtd_nodes)
return S2;
if(intercepta_esf(raio, esf_cab[id_no],&dummy)==1){
if(esf_cab[id_no].qtd_triangulos <= MIN_TRIANGULOS){
//printf("entreoi id_no=%ld\n",id_no);
for(k=0;k<esf_cab[id_no].qtd_triangulos; k++){ //Procura triângulos.
if(intercepta(raio, esf_cab[id_no].lista_triang[k] ,&S1)==1 && S1.dist<S2.dist){
S2=S1;
S2.id=id_no;
S2.id_tri=k;
}
}
return S2;
}
else{
S1=get_esf_node_simpl(raio,i_esq,esf_cab);
if(S1.dist!=TAM_RAIO+100) return S1;
S2=get_esf_node_simpl(raio,i_dir,esf_cab);
if(S2.dist!=TAM_RAIO+100) return S2;
}
}
else
return S2;
}*/
/*
* Função recursiva utilizada para fazer o trace do raio. Recebe como parâmetro o raio
* que será rastreado e um limite. O limite serve para limitar a recursão, do contrário
* o algoritmo poderia entrar em loop infinito. Primeiramente verifica se o raio intercepta algum
* triângulo presente na spheere-tree através da função get_esf_node, e depois verifica se o raio
* intercepta alguma esfera através da função intercepta_esf.
* O objeto que receberá a contribuição de cor será o objeto que estiver mais perto da origem do raio.
*/
__device__ color trace(ray raio, int lim, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes){
color col;
int k, v_intercepta=0, l;
esf aux2; triang aux1;
ipoint S1,S2;
S2.dist=TAM_RAIO+100;
if(lim >= LIMITE) {
col.r = col.g = col.b = 0.0;
return col;
}
//S2=get_esf_node(raio,0,esf_cab,qtd_nodes);
v_intercepta=(S2.dist < TAM_RAIO+100);
/*for(k=0;k<qtd_nodes;k++){
if(esf_cab[k].qtd_triangulos<=MIN_TRIANGULOS){
for(l=0;l<esf_cab[k].qtd_triangulos;l++){
if(intercepta(raio, esf_cab[k].lista_triang[l] ,&S1)==1 && S1.dist<S2.dist){
S2=S1;
aux1=esf_cab[k].lista_triang[l];
v_intercepta=1;
}
}
}
}*/
for(k=0;k<ind;k++){ //Procura esferas.
if(intercepta_esf(raio,esferas[k],&S1)==1 && S1.dist<S2.dist){
S2=S1;
aux2=esferas[k];
v_intercepta=2;
}
}
if(v_intercepta==1){
//printf("dist=%f\n",S3.dist);
col=shade(esf_cab[S2.id].lista_triang[S2.id_tri],S2,lim,raio,luz,ind,esferas,esf_cab, qtd_nodes);
}
else if(v_intercepta==2){
col=shade_esf(aux2,S2,lim,raio,luz,ind,esferas,esf_cab, qtd_nodes);
}
else {
col.r = col.g = col.b = 0.0;
}
return col;
}
/*
* Função utilizada para colorir o ponto interceptado quando este pertence a um triangluo.
* Ela recebe como parâmetro o triângulo interceptado, o ipoint que é os dados referentes ao
* ponto interceptado, o valor atual da chamada recursiva.
* A função utiliza o modelo phong para calcular a cor do pixel, ela também verifica se existe alguma
* obstrução entre o objeto interceptado e a luz para criar somra. Por fim verifica se o objeto
* em questão possui prpriedade reflexiva. Se sim, chama o trace novamente para rastrear o raio
* refletido.
* A função retorna a variável cor que serve como contribuição para um pixel da imagem.
*/
__device__ color shade(triang tri, ipoint S1, int lim, ray raio, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes){
int i, v_shadow=0, k;
obj3d luz_dir; FLOAT luz_dist;
ipoint S2; S2.dist=TAM_RAIO+100;
ray shadow_ray;
color Id, Is, cor;
Id.r = Id.g = Id.b = Is.r = Is.g = Is.b = cor.r = cor.g = cor.b = 0;
for(i=0; i<QTD_LUZ; i++){
luz_dir.x= luz[i].pos.x - S1.pos.x;
luz_dir.y= luz[i].pos.y - S1.pos.y;
luz_dir.z= luz[i].pos.z - S1.pos.z;
luz_dist=sqrt(DOT(luz_dir,luz_dir));
shadow_ray.orig = S1.pos;
shadow_ray.dir = luz_dir;
S2=get_esf_node(shadow_ray,0,esf_cab,qtd_nodes);
v_shadow=(S2.dist < luz_dist);
if(v_shadow==0){
for(k=0;k<ind;k++){
if(intercepta_esf(shadow_ray,esferas[k],&S2)==1 && S2.dist < luz_dist && S2.dist > ERR_MARGIN){
v_shadow=1;
break;
}
}
}
if(!v_shadow) {
NORMALIZE(luz_dir);
Id.r = tri.dif * (MAX(DOT(S1.normal, luz_dir), 0.0));
Is.r = tri.spec * (tri.spec > 0.0 ? pow(MAX(DOT(S1.vrefl, luz_dir), 0.0), tri.shi) : 0.0);
Id.g = Id.b = Id.r; Is.g = Is.b = Is.r;
Id.r *= luz[i].dif.r; Is.r *= luz[i].spec.r;
Id.g *= luz[i].dif.g; Is.g *= luz[i].spec.g;
Id.b *= luz[i].dif.b; Is.b *= luz[i].spec.b;
cor.r += Id.r * tri.col.r + Is.r;
cor.g += Id.g * tri.col.g + Is.g;
cor.b += Id.b * tri.col.b + Is.b;
}
if(tri.refl > 0.0) {
ray raio2;
color rcol;
raio2.orig = S1.pos;
raio2.dir = S1.vrefl;
raio2.dir.x *= TAM_RAIO;
raio2.dir.y *= TAM_RAIO;
raio2.dir.z *= TAM_RAIO;
rcol = trace(raio2, lim + 1,luz, ind, esferas, esf_cab, qtd_nodes);
cor.r += rcol.r * tri.refl;
cor.g += rcol.g * tri.refl;
cor.b += rcol.b * tri.refl;
}
}
return cor;
}
__device__ color shade_esf(esf esfera, ipoint S1, int lim, ray raio, light luz[QTD_LUZ], long ind, esf *esferas, esf *esf_cab, long qtd_nodes){
color col;
int i, v_shadow=0, k;
obj3d luz_dir; FLOAT luz_dist;
ipoint S2; S2.dist=TAM_RAIO+100;
ray shadow_ray;
color Id, Is;
Id.r = Id.g = Id.b = Is.r = Is.g = Is.b = col.r = col.g = col.b = 0;
for(i=0; i<QTD_LUZ; i++){
luz_dir.x=luz[i].pos.x - S1.pos.x;
luz_dir.y=luz[i].pos.y - S1.pos.y;
luz_dir.z=luz[i].pos.z - S1.pos.z;
luz_dist=sqrt(DOT(luz_dir,luz_dir));
shadow_ray.orig = S1.pos;
shadow_ray.dir = luz_dir;
S2=get_esf_node(shadow_ray,0,esf_cab,qtd_nodes);
v_shadow=(S2.dist < luz_dist);
if(v_shadow==0){
for(k=0;k<ind;k++){
if(intercepta_esf(shadow_ray,esferas[k],&S2)==1 && S2.dist < luz_dist && S2.dist > ERR_MARGIN){
v_shadow=1;
break;
}
}
}
if(v_shadow==0) {
NORMALIZE(luz_dir);
Id.r = esfera.dif * (MAX(DOT(S1.normal, luz_dir), 0.0));
Is.r = esfera.spec * (esfera.spec > 0.0 ? pow(MAX(DOT(S1.vrefl, luz_dir), 0.0), esfera.shi) : 0.0);
Id.g = Id.b = Id.r; Is.g = Is.b = Is.r;
Id.r *= luz[i].dif.r; Is.r *= luz[i].spec.r;
Id.g *= luz[i].dif.g; Is.g *= luz[i].spec.g;
Id.b *= luz[i].dif.b; Is.b *= luz[i].spec.b;
col.r += Id.r * esfera.col.r + Is.r;
col.g += Id.g * esfera.col.g + Is.g;
col.b += Id.b * esfera.col.b + Is.b;
//col.r=1;
//col.g=1;
//col.b=1;
}
}
if(esfera.refl > 0.0) {
ray raio2;
color rcol;
raio2.orig = S1.pos;
raio2.dir = S1.vrefl;
raio2.dir.x *= TAM_RAIO;
raio2.dir.y *= TAM_RAIO;
raio2.dir.z *= TAM_RAIO;
rcol = trace(raio2, lim + 1, luz, ind, esferas, esf_cab, qtd_nodes);
col.r += rcol.r * esfera.refl;
col.g += rcol.g * esfera.refl;
col.b += rcol.b * esfera.refl;
}
return col;
}
/*
* Função que inicia o processo de renderização. Ela recebe como parâmetro
* o valor da largura WIDTH (x) da imagem, o valor da altura HEIGH (y) da imagem
* e a quantidade de pixels que cada saíra de cada pixel. Ela chamará as funções
* trace que fará o ray tracing, e a função get_p_ray que calculará o raio primário.
* Por fim, receberá a cor final do processo, e amrazenará suas componentes na matriz
* imagem, que é uma variável global usada para escrever os valores dos pixels no
* arquivo ppm.
*/
__global__ void renderiza(int x, int y, int ray_per_pixel, long ind, esf *esferas, esf *esf_cab, color *imagem, long qtd_nodes){
long indx, i,j;
color col;
col.r=0; col.g=0; col.b=0;
light luz[QTD_LUZ];
FLOAT m[12]; ray raio;
cria_luz(0, 1,1,1, 1,1,1, 0,498,0, luz);
posiciona_camera(300,100,300, 0,0,0, 0,1,0, 1.56746300581, m);
j = (blockIdx.x * blockDim.x) + threadIdx.x;
i = (blockIdx.y * blockDim.y) + threadIdx.y;
indx=i * WIDTH + j;
__syncthreads();
if(j<WIDTH && i<HEIGH){
raio=get_p_ray(j,i,j, m);
col=trace(raio,0, luz, ind, esferas, esf_cab, qtd_nodes);
imagem[indx].r = MIN(col.r,1) * 255;
imagem[indx].g = MIN(col.g,1) * 255;
imagem[indx].b = MIN(col.b,1) * 255;
}
}
#define pivot_index() (begin+(end-begin)/2)
#define swap(a,b,t) ((t)=(a),(a)=(b),(b)=(t))
void sort(triang array[], long begin, long end, int mod) {
/*** Use of static here will reduce memory footprint, but will make it thread-unsafe ***/
static triang pivot;
static triang t; /* temporary variable for swap */
if (end > begin) {
long l = begin + 1;
long r = end;
swap(array[begin], array[pivot_index()], t); /*** choose arbitrary pivot ***/
pivot = array[begin];
while(l < r) {
switch(mod){
case 0:
if (array[l].baricentro.x <= pivot.baricentro.x) {
l++;
} else {
while(l < --r && array[r].baricentro.x >= pivot.baricentro.x) /*** skip superfluous swaps ***/
;
swap(array[l], array[r], t);
}
break;
case 1:
if (array[l].baricentro.y <= pivot.baricentro.y) {
l++;
} else {
while(l < --r && array[r].baricentro.y >= pivot.baricentro.y) /*** skip superfluous swaps ***/
;
swap(array[l], array[r], t);
}
break;
case 2:
if (array[l].baricentro.z <= pivot.baricentro.z) {
l++;
} else {
while(l < --r && array[r].baricentro.z >= pivot.baricentro.z) /*** skip superfluous swaps ***/
;
swap(array[l], array[r], t);
}
break;
}
}
l--;
swap(array[begin], array[l], t);
sort(array, begin, l, mod);
sort(array, r, end, mod);
}
}
#undef swap
#undef pivot_index
esf bound_ball(long begin, long end, long qtd_triangulos){
obj3d C, dPx, dPy, dPz, dP, Paux;
FLOAT rad, rad2, dx2, dy2, dz2, dist, dist2;
FLOAT xmin, xmax, ymin, ymax, zmin, zmax;
long Pxmin, Pxmax, Pymin, Pymax, Pzmin, Pzmax;
long qtd_pontos = qtd_triangulos * 3;
obj3d *P;
long i,j=0;
esf B;
P= (obj3d*) malloc(sizeof(obj3d) * qtd_pontos);
for(i=begin;i<=end;i++){
P[j]=faces[i].A;
P[j+1]=faces[i].B;
P[j+2]=faces[i].C;
j+=3;
}
xmin = xmax = P[0].x;
ymin = ymax = P[0].y;
zmin = zmax = P[0].z;
Pxmin = Pxmax = Pymin = Pymax = Pzmin = Pzmax = 0;
for(i=0;i<j;i++){
if (P[i].x < xmin) {
xmin = P[i].x;
Pxmin = i;
}
else if (P[i].x > xmax) {
xmax = P[i].x;
Pxmax = i;
}
if (P[i].y < ymin) {
ymin = P[i].y;
Pymin = i;
}
else if (P[i].y > ymax) {
ymax = P[i].y;
Pymax = i;
}
if (P[i].z < zmin) {
zmin = P[i].z;
Pzmin = i;
}
else if (P[i].z > zmax) {
zmax = P[i].z;
Pzmax = i;
}
}
dPx.x = P[Pxmax].x - P[Pxmin].x;
dPx.y = P[Pxmax].y - P[Pxmin].y;
dPx.z = P[Pxmax].z - P[Pxmin].z;
dPy.x = P[Pymax].x - P[Pymin].x;
dPy.y = P[Pymax].y - P[Pymin].y;
dPy.z = P[Pymax].z - P[Pymin].z;
dPz.x = P[Pzmax].x - P[Pzmin].x;
dPz.y = P[Pzmax].y - P[Pzmin].y;
dPz.z = P[Pzmax].z - P[Pzmin].z;
dx2 = DOT(dPx,dPx);
dy2 = DOT(dPy,dPy);
dz2 = DOT(dPz,dPz);
if (dx2 >= dy2 && dx2 >= dz2) {
C.x = P[Pxmin].x + (dPx.x / 2.0);
C.y = P[Pxmin].y + (dPx.y / 2.0);
C.z = P[Pxmin].z + (dPx.z / 2.0);
Paux.x=P[Pxmax].x - C.x;
Paux.y=P[Pxmax].y - C.y;
Paux.z=P[Pxmax].z - C.z;
rad2 = DOT(Paux,Paux);
}
else if (dy2 >= dx2 && dy2 >= dz2){
C.x = P[Pymin].x + (dPy.x / 2.0);
C.y = P[Pymin].y + (dPy.y / 2.0);
C.z = P[Pymin].z + (dPy.z / 2.0);
Paux.x=P[Pymax].x - C.x;
Paux.y=P[Pymax].y - C.y;
Paux.z=P[Pymax].z - C.z;
rad2 = DOT(Paux,Paux);
}
else {
C.x = P[Pzmin].x + (dPz.x / 2.0);
C.y = P[Pzmin].y + (dPz.y / 2.0);
C.z = P[Pzmin].z + (dPz.z / 2.0);
Paux.x=P[Pzmax].x - C.x;
Paux.y=P[Pzmax].y - C.y;
Paux.z=P[Pzmax].z - C.z;
rad2 = DOT(Paux,Paux);
}
rad = sqrt(rad2);
for (i=0; i<j; i++) {
dP.x = P[i].x - C.x;
dP.y = P[i].y - C.y;
dP.z = P[i].z - C.z;
dist2 = DOT(dP, dP);
if (dist2 <= rad2)
continue;
dist = sqrt(dist2);
rad = (rad + dist) / 2.0;
rad2 = rad * rad;
C.x = C.x + ((dist-rad)/dist) * dP.x;
C.y = C.y + ((dist-rad)/dist) * dP.y;
C.z = C.z + ((dist-rad)/dist) * dP.z;
}
B=construct_esf(C.x,C.y,C.z,rad, 0,1,0,1, 1,1,1,0);
B.qtd_triangulos=qtd_triangulos;
free(P);
return B;
}
int cria_arv_esf(long ini, long fim, long qtd_triang, int mod, long id_no){
int i,j=0;
long aux=fim - ini + 1;
/*if(ini>=fim){
printf("DEU RUIM\n\n\n\n\n");
return 1;
}*/
if(qtd_triang<=MIN_TRIANGULOS){
h_esf_cab[id_no]=bound_ball(ini,fim, aux);
for(i=ini;i<=fim;i++){
h_esf_cab[id_no].lista_triang[j]=faces[i];
j++;
}
return 0;
}
sort(faces,ini,fim,mod);
//printf("sorted...\n");
h_esf_cab[id_no]=bound_ball(ini,fim,aux);
//printf("balls\n");
mod = mod % 2;
mod++;
qtd_triang=aux/2;
cria_arv_esf(ini, fim - qtd_triang, qtd_triang, mod, 2 * (id_no + 1) - 1);
cria_arv_esf(ini + qtd_triang,fim, qtd_triang, mod, 2 * (id_no + 1));
return 0;
}
void montar_chao(long *ntriangles){
obj3d A1,B1,C1;
long paredes=6;
long triangulos=paredes*2;
FLOAT r_parede=0.5, g_parede=0.4, b_parede=0.6, a_parede=0;
FLOAT spec_parede=0.0, dif_parede=1, shin_parede=0.0, refl_parede=1;
FLOAT r_chao=1, g_chao=1, b_chao=1, a_chao=0;
FLOAT spec_chao=0, dif_chao=1, shin_chao=0, refl_chao=0;
FLOAT r_teto=1, g_teto=0, b_teto=0, a_teto=0;
FLOAT spec_teto=0, dif_teto=1, shin_teto=0, refl_teto=0;
*ntriangles+=triangulos ;
faces = (triang*) malloc(sizeof(triang) * (*ntriangles));
jind=triangulos;
//chão
A1.x=-1000; A1.y=-2; A1.z=-1000;
B1.x=-1000; B1.y=-2; B1.z=1000;
C1.x=1000; C1.y=-2; C1.z=-1000;
faces[0]=construc_triang(A1,B1,C1, r_chao, g_chao, b_chao, a_chao, spec_chao, dif_chao, shin_chao, refl_chao);
faces[1]=construc_triang(escala(A1,-1,1,-1),escala(B1,-1,1,-1),escala(C1,-1,1,-1),
r_chao, g_chao, b_chao, a_chao, spec_chao, dif_chao, shin_chao, refl_chao);
//parede 2
A1.x=-1000; A1.y=-1000; A1.z=500;
B1.x=1000; B1.y=1000; B1.z=500;
C1.x=1000; C1.y=-1000; C1.z=500;
faces[2]=construc_triang(A1,B1,C1, r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
faces[3]=construc_triang(escala(A1,-1,-1,1),escala(B1,-1,-1,1),escala(C1,-1,-1,1),
r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
//parede 3
A1.x=-500; A1.y=-1000; A1.z=-1000;
B1.x=-500; B1.y=1000; B1.z=-1000;
C1.x=-500; C1.y=1000; C1.z=1000;
faces[4]=construc_triang(A1,B1,C1, r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
faces[5]=construc_triang(escala(A1,1,-1,-1),escala(B1,1,-1,-1),escala(C1,1,-1,-1),
r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
//teto
A1.x=-1000; A1.y=500; A1.z=-1000;
B1.x=-1000; B1.y=500; B1.z=1000;
C1.x=1000; C1.y=500; C1.z=-1000;
faces[6]=construc_triang(B1,A1,C1, r_teto, g_teto, b_teto, a_teto, spec_teto, dif_teto, shin_teto, refl_teto);
faces[7]=construc_triang(escala(B1,-1,1,-1),escala(A1,-1,1,-1),escala(C1,-1,1,-1),
r_teto, g_teto, b_teto, a_teto, spec_teto, dif_teto, shin_teto, refl_teto);
//parede 2 - complemento
A1.x=-1000; A1.y=-1000; A1.z=-500;
B1.x=1000; B1.y=1000; B1.z=-500;
C1.x=1000; C1.y=-1000; C1.z=-500;
faces[8]=construc_triang(B1,A1,C1, r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
faces[9]=construc_triang(escala(B1,-1,-1,1),escala(A1,-1,-1,1),escala(C1,-1,-1,1),
r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
//parede 3 - complemento
A1.x=500; A1.y=-1000; A1.z=-1000;
B1.x=500; B1.y=1000; B1.z=-1000;
C1.x=500; C1.y=1000; C1.z=1000;
faces[10]=construc_triang(B1,A1,C1, r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
faces[11]=construc_triang(escala(B1,1,-1,-1),escala(A1,1,-1,-1),escala(C1,1,-1,-1),
r_parede, g_parede, b_parede, a_parede, spec_parede, dif_parede, shin_parede, refl_parede);
}
static int vertex_cb(p_ply_argument argument) {
long eol;
ply_get_argument_user_data(argument, NULL, &eol);
switch (eol) {
case 0:
vertices[jind_v].x=ply_get_argument_value(argument);
break;
case 1:
vertices[jind_v].y=ply_get_argument_value(argument);
break;
case 2:
vertices[jind_v].z=ply_get_argument_value(argument);
jind_v++;
break;
default:
break;
}
return 1;
}
static int face_cb(p_ply_argument argument) {
long length, value_index;
ply_get_argument_property(argument, NULL, &length, &value_index);
switch (value_index) {
case 0:
vert_aux1=ply_get_argument_value(argument);
//printf("caso %ld Indice %ld\n", value_index, vert_aux1);
break;
case 1:
vert_aux2=ply_get_argument_value(argument);
//printf("caso %ld Indice %ld\n", value_index, vert_aux2);
break;
case 2:
vert_aux3=ply_get_argument_value(argument);
faces[jind]=construc_triang(trans(escala(vertices[vert_aux1],1000,1000,1000),0,30,0),
trans(escala(vertices[vert_aux2],1000,1000,1000),0,30,0),
trans(escala(vertices[vert_aux3],1000,1000,1000),0,30,0),
1,1,1,1 ,0,1,1,0);
jind++;
break;
default:
break;
}
return 1;
}
int ler_ply(char *fname) {
long nvertices, ntriangles;
jind_v=0; jind=0;
p_ply ply = ply_open(fname, NULL, 0, NULL);
if (!ply) return 1;
if (!ply_read_header(ply)) return 1;
nvertices = ply_set_read_cb(ply, "vertex", "x", vertex_cb, NULL, 0);
ply_set_read_cb(ply, "vertex", "y", vertex_cb, NULL, 1);
ply_set_read_cb(ply, "vertex", "z", vertex_cb, NULL, 2);
ntriangles = ply_set_read_cb(ply, "face", "vertex_indices", face_cb, NULL, 0);
vertices = (obj3d*) malloc(sizeof(obj3d) * nvertices);
montar_chao(&ntriangles);
printf("%ld\n%ld\n", nvertices, ntriangles);
if (!ply_read(ply)) return 1;
ply_close(ply);
free(vertices);
jind--;
return 0;
}
long ler_arq_esferas(){
int i;
long ind=1;
h_esferas=(esf*) malloc(sizeof(esf) * ind);
for(i=0;i<ind;i++){
h_esferas[i]=construct_esf(200,100,0,50, 0,0,1,1 ,0,1,0,0);
}
return ind;
}
long get_qtd_nos(){
long i=jind, h=0, qtd_nos;
while(i>=10){
i/=2;
h++;
}
altura=h;
qtd_nos= (pow(2,h) - 1) * 2 +1;
h_esf_cab = (esf*) calloc(qtd_nos,sizeof(esf));
printf("\nqtd_nos=%ld h=%ld\n", qtd_nos, h);
return qtd_nos;
}
int main(int argc, char** argv){
long ind, qtd_nos;
h_imagem=(color*) calloc(HEIGH * WIDTH,sizeof(color));
ler_ply("bun_zipper.ply");
ind=ler_arq_esferas();
qtd_nos=get_qtd_nos();
cria_arv_esf(0,jind,jind,0,0);
dim3 threadsPerBlock(16, 16);
dim3 numBlocks((WIDTH/threadsPerBlock.x)+1, (HEIGH/threadsPerBlock.y)+1);
CUDA_CHECK_RETURN(cudaMalloc((void**) &d_esf_cab, qtd_nos * sizeof(esf)));
CUDA_CHECK_RETURN(cudaMalloc((void**) &d_esferas, ind * sizeof(esf)));
CUDA_CHECK_RETURN(cudaMalloc((void**) &imagem, sizeof(color) * HEIGH * WIDTH));
CUDA_CHECK_RETURN(cudaMemcpy(d_esf_cab, h_esf_cab, qtd_nos * sizeof(esf), cudaMemcpyHostToDevice));
CUDA_CHECK_RETURN(cudaMemcpy(d_esferas, h_esferas, ind * sizeof(esf), cudaMemcpyHostToDevice));
printf("comecei a renderizar\n");
renderiza <<<numBlocks,threadsPerBlock>>> (WIDTH,HEIGH,1, ind, d_esferas, d_esf_cab, imagem, qtd_nos);
printf("terminei de renderizar\n");
cudaDeviceSynchronize();
CUDA_CHECK_RETURN(cudaMemcpy(h_imagem, imagem, HEIGH * WIDTH * sizeof(color), cudaMemcpyDeviceToHost));
printf("comecei a criar imagem\n");
create_image_ppm("teste13_64_cuda.ppm");
printf("terminei criar imagem\n");
cudaFree(d_esferas);
cudaFree(d_esf_cab);
cudaFree(imagem);
free(h_imagem);
free(h_esf_cab);
free(faces);
free(h_esferas);
//sort(faces,0,jind,2);
//bound_ball(0,jind, jind, esf_cab);
/*for(i=0;i<qtd_nodes;i++){
printf("esf_node[%d].qtd_triangulos=%ld raio=%f centro(%f,%f,%f)\n",i,esf_cab[i].qtd_triangulos,esf_cab[i].r,
esf_cab[i].centro.x, esf_cab[i].centro.y, esf_cab[i].centro.z);
}*/
/*for(i=0; i<jind; i++)
dump_triang(faces[i],i);*/
//printf("esfera: centro(%f, %f, %f) raio:%f\n", esf_cab->centro.x, esf_cab->centro.y, esf_cab->centro.z, esf_cab->r);
return 0;
}
On my main function I do the allocations for the GPU, can you see something wrong there?
Thank you very much for the help you have be given to me.