Hey guys, could you help me out. I don’t know what is wrong with my code, I tried every thing, but the error of illegal memory acces continues.
If it is not to much trouble to ask, can you help me find the error on this one?
This is my code:
/* *
* 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 3
#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;
long 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;
}
/*
* 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;
}
else{
//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;
S2.id=k;
S2.id_tri=l;
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 ray_per_pixel, long ind, esf *esferas, esf *esf_cab, color *imagem, long qtd_nodes){
long indx, i,j;
color col;
col.r=0.0; col.g=0.0; col.b=0.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;
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=0;
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=2;
h_esferas=(esf*) malloc(sizeof(esf) * ind);
//for(i=0;i<ind;i++){
h_esferas[0]=construct_esf(200,100,0,50, 0,0,1,1 ,0,1,0,0);
h_esferas[1]=construct_esf(0,0,0,200, 0,0,1,1 ,1,1,1,1);
//}
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), (HEIGH/threadsPerBlock.y));
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>>> (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);
return 0;
}
Thanks for the help. Good night to everyone.