/* Last edited on 2008-08-01 03:09:39 by stolfi */ //rdo_synthesis.c { /* Compute the direction of the camera ray through the image center: */ vector3_t ctr_dir; rdo_camera_compute_ray_direction(0.0, 0.0, cam, &ctr_dir); /* Get the image pixel site that is nearest to the image center: */ int ctr_pix_index = rdo_raytrace_find_nearest_visible(sds, &(cam->eye), &ctr_dir, pix); site_t *ctr_pix_s = (ctr_pix_index < 0 ? NULL : &(pix->e[ctr_pix_index])); fprintf(stderr, "central pixel site = "); if (ctr_pix_s == NULL) { fprintf(stderr, "NULL"); } else { rdo_site_write(stderr, ctr_pix_s); } fprintf(stderr, "\n"); } void rdo_synthesis_compute_coeffs_realistic_image ( Opcoes *o, camera_t *cam, scene_t *scn, image_t *img, Luz luzes[], int num_luzes, approx_basis_t *bas, int channel, int quicadas, double L[] ) { // !!! Simplificar -- transformar luzes em objetos com emissão própria. !!! //double pre_ini, pre_fim; //instantes de início e fim do pré-processamento. //double mul_ini, mul_fim; //instantes de início e fim da iteração de rdo_radiosity. //calcula matriz {D} de transf. de radiância de luzes externas para pontos de amostra: dspmat_t *D = ObtemMatrizDeTransferenciaDasLuzes(o, scn, channel, bas); //calcula radiâncias amostradas {B[i]} devidas a emissão própria e luzes externas: double *B = NOTNULL(malloc(sizeof(double)*bas->nel)); CalculaRadianciasIniciais(luzes, num_luzes, scn, channel, bas, &D, B); //converte as radiâncias amostradas {B[i]} para coeficientes {E[j]} da base: double *E = NOTNULL(malloc(sizeof(double)*bas->nel)); ConverteRadianciasAmostradasParaCoefsDaBase(&(bas->MM_Inv), B, E, bas->nel); //inicializa os coefs {L} da radiância total com os coefs iniciais {E}: InicializaRadianciaDaBase(E, L, bas->nel); if(quicadas > 0) { //obtém a matriz {T} de transferência de radiância entre os sites de amostragem: dspmat_t T; ObtemMatrizDeTransferenciaDosCentros(o, scn, channel, bas, &T); //constrói a matriz {R} de transferência dos coeficientes da base: dspmat_t R; CalculaMatrizDeTransferenciaDaBase(scn, channel, bas, &T, &(bas->MM), &(bas->MM_Inv), &R); //iteração da radiosidade: // time(&mul_ini); int iq; for(iq = 1; iq <= o->quicadas; iq++) { fprintf(stderr, "cálculo de radiosidade -- iteração %d\n", iq); IteracaoDeTransferenciaDaBase(&R, E, L, bas->nel); } // time(&mul_fim); Desalocadspmat_t(&R); Desalocadspmat_t(&T); } free(B); free(E); Desalocadspmat_t(&D); } //rdo_scene.c dspmat_t rdo_scene_get_site_transfer_matrix(scene_t *scn, int channel) { //obtem a matriz de transferência {scn->T} para albedo unitário: if (scn->T.ents == 0) { scn->T = rdo_lighting_compute_site_transfer_matrix ( &(scn->sds), MAX_ALBEDO, scn->smp.e, scn->smp.ne, scn->smp.e, scn->smp.ne ); } //combina com os albedos do channel {c}. dspmat_t T = rdo_lighting_combine_transfer_matrix_and_diffusion_coeffs (&(scn->T), MAX_ALBEDO, &(scn->sds), &(scn->fns), channel); return T; } //rdo_linear_algebra.c if (Aii == 0.0) { fprintf(stderr, "%s:%d: ** zero element in diagonal\n", __FILE__, __LINE__); fprintf(stderr, " row [%d] = (", i); uint32_t p; for (p = 0; p < A->ents; p++) { dspmat_entry_t *Ap = &(A->e[p]); if (Ap->row == i) { fprintf(stderr, " [%d]=%24.16e", Ap->col, Ap->val); } } fprintf(stderr, " )\n"); } //rdo_synthesis.h dspmat_t rdo_synthesis_compute_rendering_matrix(approx_basis_t *bas, site_vec_t *pix); /* Computes the matrix {Q} that converts a column vector of coefficients of the basis {bas} to the radiances of the pixel sites in the list {pix}. */ dspmat_t rdo_synthesis_compute_rendering_matrix(approx_basis_t *bas, site_vec_t *pix) { int n_bas = bas->ne; int n_pix = pix->ne; dspmat_t Q = dspmat_new(n_pix, n_bas, 10*(n_bas + n_pix)); double val[n_bas]; /* Values of all basis elements at a given pixel. */ uint32_t posQ = 0; int ip; for (ip = 0; ip < n_pix; ip++) { site_t *sP = &(pix->e[ip]); rdo_approx_basis_eval_all_elements(bas, sP, val); posQ = dspmat_add_row(&Q, posQ, ip, val, n_bas); } dspmat_trim(&Q, posQ); return Q; } void rdo_synthesis_complete_scene ( scene_t *scn, camera_t *cam, dist_type_t tpDist, basis_elem_type_t tpElem, sampling_options_t *opSmp, approx_basis_options_t *opBas ); /* Creates any auxiliary component of the scene that is missing or empty. If the sampling sites {scn->smp} are missing, recreates them with the {opSmp} options. If the approx basis {scn->bas} is missing, recreates it from the sampling sites with the {opBas} options. If the {T} matrix is missing, recomputes it from the {scn->smp} sites. If the {M} matrix is missing, recomputes it from the {scn->bas} basis. */ //rdo_image.h //um pixel do img fotografico: typedef struct Pixel { point3_t center; frgb_t cor; } Pixel; image_t *rdo_image_new(int channels, int nh, int nv); /* Cria um img virtual (imagem) com {nh} colunas e {tamanho.c[1]} linhas. */ void rdo_image_free(image_t *img); void rdo_image_set_sample(image_t *img, int channel, int ih, int iv, float valor); /* Coloca {valor} na amostra {channel} do pixel de índices {ih,iv} no img {img}. */ void rdo_image_set_pixel(image_t *img, int ih, int iv, frgb_t *cor); /* Pinta o pixel de índices {ih,iv} no img {img} com a {cor} dada. */ image_t *rdo_image_new(int channels, int nh, int nv) { return float_image_new(channels, nh, nv); } void rdo_image_free(image_t *img) { float_image_free(img); /* Discards the header, too. */ } void rdo_image_draw_cross(image_t *img, int ih, int iv, frgb_t *cor, int radius, bool_t diagonal) { rdo_image_set_pixel(img, ih, iv, cor); int k; for (k = 1; k <= radius; k++) { if (diagonal) { rdo_image_set_pixel(img, ih+k, iv+k, cor); rdo_image_set_pixel(img, ih+k, iv-k, cor); rdo_image_set_pixel(img, ih-k, iv+k, cor); rdo_image_set_pixel(img, ih-k, iv-k, cor); } else { rdo_image_set_pixel(img, ih+k, iv, cor); rdo_image_set_pixel(img, ih-k, iv, cor); rdo_image_set_pixel(img, ih, iv+k, cor); rdo_image_set_pixel(img, ih, iv-k, cor); } } } void rdo_image_set_pixel(image_t *img, int ih, int iv, frgb_t *cor) { int nc = img->sz[0]; int nh = img->sz[1]; int nv = img->sz[2]; int channel; for (channel = 0; channel < nc; channel++) { if ((ih >= 0) && (ih < nh) && (iv >= 0) && (iv < nv)) { float val = (channel >= 3 ? frgb_Y(cor) : cor->c[channel]); float_image_set_sample(img, channel, ih, iv, val); } } } //srdo_sampling.c void rdo_sampling_generate_sites_by_projection ( solid_vec_t *sds, camera_t *cam, int nh, int nv, site_vec_t *sts ) { //calcula espaçamento médio entre os radiuss no plano da imagem: double nr = hypot((double)nh, (double)nv); double passo = cam->radius/nr; //percorre a grade de pontos no plano da imagem: int nold = sts->ne; //número original de sites em {sts} int ntot = ntot; //numero total de sites em {sts}. int kh, kv; for(kh = 0; kh < nh; kh++) for(kv = 0; kv < nv; kv++) { //escolhe a posição {(h,v)} do radius no sistema de coordenadas da câmera: double h = (kh - ((double)nh)/2 + 0.1 + 0.8*drandom())*passo; double v = (((double)nv)/2 - kv + 0.1 + 0.8*drandom())*passo; //calcula a intersecção {p} do radius com o plano da imagem, no espaço: point3_t p = cam->ctr; r3_mix_in(+h, &cam->r, &p); r3_mix_in(-v, &cam->s, &p); //determina a direcao do radius: vector3_t dir = rdo_geometry_pt_pt_dir(&(cam->eye), &p); //traça o radius e armazena o resultado: site_t sR; double dR; rdo_raytrace_first_hit_site(sds, &(cam->eye), &dir, &sR, &dR); if (isfinite(dR) && (sR.isd >= 0)) { site_vec_expand(sts, ntot); sts->e[ntot] = sR; ntot++; } } fprintf(stderr, "gerados %d sites por projecao (total %d)\n ", ntot - nold, ntot); site_vec_trim(sts, ntot); } //rdo_synthesis.c void rdo_syntesis_compute_coeffs_simple_radiance(scene_t *scn, int channel, approx_basis_t *bas, double L[]) { int i; double *B = NOTNULL(malloc(sizeof(double)*bas->nel)); for(i = 0; i < bas->nel; i++) { site_t *si = &(bas->el[i].center); frgb_t cor = rdo_syntesis_simple_radiance_function(si); B[i] = cor.c[channel]; } Multiplicadspmat_tVetor(&(bas->MM_Inv), B, L, bas->nel); free(B); } void rdo_syntesis_compute_coeffs_show_interpolators(camera_t *cam, scene_t *scn, image_t *img, int channel, approx_basis_t *bas, double L[]) { int i; double *B = NOTNULL(malloc(sizeof(double)*bas->nel)); int INDICE1 = EncontraCentroProximo(100,100,cam,img,scn,bas); int INDICE2 = EncontraCentroProximo(190,230,cam,img,scn,bas); int INDICE3 = EncontraCentroProximo(130,140,cam,img,scn,bas); int INDICE4 = EncontraCentroProximo(200,180,cam,img,scn,bas); for(i = 0; i < bas->nel; i++) { if(i != INDICE1 && i != INDICE2 && i != INDICE2 && i != INDICE3 && i != INDICE4) { B[i] = 0.0; } else { B[i] = 1.0; } } Multiplicadspmat_tVetor(&(bas->MM_Inv), B, L, bas->nel); free(B); } void rdo_syntesis_compute_coeffs_show_basis_elems(camera_t *cam, scene_t *scn, image_t *img, int channel, approx_basis_t *bas, double L[]) { int i; int ind1 = EncontraCentroProximo(100,100,cam,img,scn,bas); int ind2 = EncontraCentroProximo(190,230,cam,img,scn,bas); int ind3 = EncontraCentroProximo(130,140,cam,img,scn,bas); int ind4 = EncontraCentroProximo(200,180,cam,img,scn,bas); for(i = 0; i < bas->nel; i++) { if(i != ind1 && i != ind2 && i != ind2 && i != ind3 && i != ind4) { L[i] = 0.0; } else { L[i] = 1.0; } } } void rdo_syntesis_compute_coeffs_simple_radiance(solid_vec_t *sds, int channel, approx_basis_t *bas, double L[]); /* Calcula os coeficientes {L[i]} na base {bas} que interpolam o channel {channel} da função de radiância {rdo_syntesis_simple_radiance_function}. */ void rdo_syntesis_compute_coeffs_show_interpolators(camera_t *cam, solid_vec_t *sds, image_t *img, int channel, approx_basis_t *bas, double L[]); /* Calcula os coeficientes da radiancia para ilustrar elementos interpoladores da base. Para isso escolhe alguns centróides interssantes, e determina os coeficientes que fazem a radiancia 1 nesses centers, 0 nos demais centers. */ void rdo_syntesis_compute_coeffs_show_basis_elems(camera_t *cam, solid_vec_t *sds, image_t *img, int channel, approx_basis_t *bas, double L[]); /* Calcula os coeficientes da radiancia para ilustrar elementos da base. Para isso escolhe alguns centróides interssantes, e define o valor dos coeficientes como sendo 1 nesses centers, 0 nos demais centers. */ void rdo_syntesis_compute_coeffs_realistic_image(camera_t *cam, scene_t *scn, image_t *img, int channel, int quicadas, double L[]); /* Calcula os coeficientes da radiancia para uma imagem realista da scene, para o channel {channel}. Usa os parametros de iluminação, numero de saltos, etc.. */ void rdo_syntesis_compute_coeffs_show_distance(camera_t *cam, image_t *img, scene_t *scn, site_t s[]); //!!! LATER: void filmeParaImagem(image_t *img, float_image_t *Img); //!!! LATER: void ex_write_image(FILE *wr, char *name, float_image_t *img); //!!! LATER: void corrigeGama(float_image_t *A, double gama); void PintaCentrosDosbasis_elem_ts(approx_basis_t *bas, scene_t *scn, camera_t *cam, image_t *img, frgb_t *cor, int radius, bool_t diagonal) { int i; for(i = 0; i < bas->nel; i++) { //pinta o centróide do elemento {i} se for visível: site_t *sti = &(bas->el[i].center); //centróide do elemento {i} da base. double vf = rdo_site_pt_site_visibility(cam->eye, sti->position, sti->normal, scn); if (vf > 0.0) { //centróide do elemento é visível: int ih, iv; rdo_camera_project_point(cam, bas->el[i].center.position, &ih, &iv); rdo_image_draw_cross(img, ih, iv, cor, radius, diagonal); } } } void rdo_syntesis_compute_coeffs_show_distance ( camera_t *cam, image_t *img, scene_t *scn, site_tsDePixeis *vis ) { //determina o site visível {sR} que está no center da imagem: vector3_t ray; ray.x = cam->scene.x - cam->eye.x; ray.y = cam->scene.y - cam->eye.y; ray.z = cam->scene.z - cam->eye.z; NormalizaVetor(&ray); frgb_t cor; site_t sR = EncontraInterseccao(scn,cam->eye,ray,&cor); if(sR.ind > 0) { //calcula radiâncias dos pixeis da imagem em função da distância a {sR}: CoeficientesDistancia(scn,sR,img,vis); } else { fprintf(stderr, "%s: não conseguiu escolher o site de referência\n", __FUNCTION__); exit(1); } } void rdo_syntesis_create_simple_radiance_image ( camera_t *cam, image_t *img, solid_vec_t *sds, approx_basis_t *bas, site_t s[] ); /* Armazena na película {img} uma imagem da scene produzida apenas pela {rdo_syntesis_simple_radiance_function}, sem levar em conta as luzes externas nem a rdo_radiosity. Também coloca em {Vis} todos os sites visíveis da scene correspondentes aos centers dos pixels da imagem. */ void rdo_syntesis_create_simple_radiance_image ( camera_t *cam, image_t *img, solid_vec_t *sds, approx_basis_t *bas, site_t s[] ) { int nh = img->sz[1]; int nv = img->sz[2]; for(ih = 0; ih < nh; ih++) for(iv = 0; iv < nv; iv++) { Pixel *p = &(img->pixel[ih+iv*img->n_colunas]); double h = ((double)ih) + 0.5; double v = ((double)iv) + 0.5; vector3_t dir; rdo_camera_compute_ray_direction(cam, h, v, nh, nv, &dir); site_t sp; /* Site visible at center of pixel. */ double dp; /* Its distance from the camera. */ rdo_raytrace_first_hit_site(sds, &cam->eye, &dir, &sp, &dp); if(s.isd > 0) { //intersecção valida, armazena em {Vis}: vis->s_i[cont] = sint; vis->h_pixel[cont] = ih; vis->v_pixel[cont] = iv; frgb_t cor = rdo_synthesis_simple_radiance_function(&s); cont++; } } vis->nvs = cont; printf("Gerados %d sites de pixeis visíveis\n",cont); } //cenario.c MatrizEsparsa *ObtemMatrizDeTransferenciaDasLuzes(Cenario *cen, int canal) { if (bas->Db == NULL) { bas->Db = NOTNULL(malloc(sizeof(MatrizEsparsa))); *(bas->Db) = CalculaMatrizDeTransferenciaDasLuzes ( cen, bas->sb, bas->nel ); } return bas->Db; } //sintese.h void DescartaItensDaBase(bool_t sitios, bool_t elems, bool_t matrizes, BaseAprox **basP); /* Se {matrizes} for {TRUE}, descarta apenas as matrizes de transferência contidas em {**basP}, se as houver. Se {elems} for {TRUE}, descarta os elementos da base de {**basP}, se os houver (e portanto também as matrizes}. Se {sitios} for {TRUE}, descarta o conjunto de sítios de amostragem de {*basP} (e portanto as matrizes e os elementos da base). Neste último caso, tamb;em descarta {*basP} e faz {*basP = NULL}. */ void DescartaItensIncompativeisDaBase ( Cena *cena, TipoDeDistancia tpDist, TipoDeElemento tpElem, OpcoesDeAmostragem *opAmo, OpcoesDeBase *opBas, BaseAprox **basP ); /* Verifica se as matrizes, elementos da base, e sítios de amostragem de {**basP} são compatíveis com as opções do usuário {tpDist}, {opAmo} e {opBas}. Itens que forem considerados incompatíveis são descartados. Nota: algumas incompatibilidades podem não ser detectadas. */ //iluminacao.h MatrizEsparsa CalculaMatrizDeTransferenciaDasLuzes(ListaDeObjetos *obs, Sitio sd[], int nsd); /* Calcula a matriz {T} de transferência de radiosidade das luzes externas da {cena} para os sítios {sd[0..nsd-1]}, para o canal {canal}. O produto desta matriz pelo vetor coluna {S[j]} de intensidades das luzes externas produz o vetor coluna {B[i]} das radiâncias nos sítios {sd[i]}, calculadas por simples ray-tracing (uma única quicada), sem levar em conta a emissão própria dos elementos da cena nem iluminação indireta. !!! Deveria receber a matriz de visibilidade já calculada. !!! Será que precisa da métrica e/ou dos raios? */ double CalculaCoefDeIluminacaoDireta(Cena *cena, Ponto3d *pL, double dL, Ponto3d *pD, Vetor3d *nD); /* Calcula a contribuição de uma luz pontual situada no ponto {pL} para a radiância da superfície da cena no ponto {pD}, cuja normal é {nD}. Considera apenas fótons dessa fonte de luz que são captados pelo olho após uma única quicada na cena. Supõe que o albedo da superfície no ponto {pD} é 1.0, e que a intensidade da luz é tal que resultaria em uma radiância unitária para uma superfície lambertiana de albedo 1 iluminada normalmente a uma distância {dL} de {pL}. */ MatrizEsparsa CalculaMatrizDeTransferenciaDasLuzes(ListaDeObjetos *obs, Sitio sd[], int nsd) { int nlz = cena->num_luzes; MatrizEsparsa D = AlocaMatrizEsparsa(nsd*nlz, nsd, nlz); int i; for(i = 0; i < nsd; i++) { //pega posição e normal do sítio de destino número {i}: Ponto3d *pdi = &(sd[i].posicao); Vetor3d *ndi = &(sd[i].normal); int j; for(j = 0; j < nlz; j++) { //pega as coordenadas da luz externa número {j}: Ponto3d *plj = &(cena->luzes[j].posicao); double dlj = cena->luzes[j].d_ref; //calcula e armazena o elemento {D[i,j]} da matriz. double Dij = CalculaCoefDeIluminacaoDireta(cena, plj, dlj, pdi, ndi); if (Dij != 0.0) { InsereElementoEmMatrizEsparsa(&D, Dij, i, j); } } } return D; } double CalculaCoefDeIluminacaoDireta(Cena *cena, Ponto3d *pL, double dL, Ponto3d *pD, Vetor3d *nD) { double vis = VisibilidadePontoParaSitio(cena, pL, pD, nD); //0 se invisível, 1 se visível. if (vis > 0.0) { //leva em conta a distância e inclinação de {sB}: Vetor3d rSL = DirecaoDePontoParaPonto(pD, pL); double cSL = r3_dot(nD, &rSL); if (cSL <= 0.0) { return 0.0; } double dSL = r3_dist(pL, pD) / dL; return vis*cSL/(dL*dL); } else { return 0.0; } } double CalculaCoefRadianciaEntreSitios(Cena *cena, Ponto3d *pO, Vetor3d *nO, Ponto3d *pD, Vetor3d *nD) { double vis = VisibilidadeSitioParaSitio(cena, pO, nO, pD, nD); //0 se invisível, 1 se visível. if (vis > 0.0) { double valor = 0.0; //pertubacao das coordenadas do centro i pn.x += ALPHA*nj.x; pn.y += ALPHA*nj.y; pn.z += ALPHA*nj.z; //vetor na direcao luz p/ centro i ray.x = pn.x - luz.x; ray.y = pn.y - luz.y; ray.z = pn.z - luz.z; NormalizaVetor(&ray); s = EncontraInterseccao(cena,luz,ray,&c,&indice,&tipo); if(s > 0.0) { p = RetornaPontonInterseccao(luz,ray,s); //se o primeiro ponto interceptado é o centro j if(Distancia2Pontos(p,bas->el[j].centro) < ERRO) { //encontra a normal do centro j if(tipo == BOLA) EncontraNormalBola(cena->bolas[indice],bas->el[j].centro,&nj); else EncontraNormalCubo1(cubos[indice],bas->el[j].centro,&nj,bas,Vis,Cam.l[j]); NormalizaVetor(&nj); //encontra a normal do centro i if(bas->sb[i].tipo == BOLA) EncontraNormalBola(cena->bolas[bas,Vis,Cam.ind[i]],bas->sb[i],&normal); else EncontraNormalCubo1(cubos[bas,Vis,Cam.ind[i]],bas->sb[i],&normal,bas,Vis,Cam.l[i]); NormalizaVetor(&normal); //verifica se o fator geometrico é valido fat1 = ProdutoEscalar(ray,normal); if(fat1 < 0.0) { fat1 = 0.0; } InverteVetor(&ray); fat2 = ProdutoEscalar(ray,nj); if(fat2 < 0.0) { fat2 = 0.0; } dist = Distancia2Pontos(bas->sb[i],bas->el[j].centro); raio = (1.0)/sqrt(Alpha[j]); area = (PI*pow(raio,2))/12; fat3 = (area/(PI*pow(dist,2))); //if(dist < 1.0) { printf("ERRO TOTAL"); scanf("%d",&i);} fator = fat1*fat2*fat3; if(i == j) fator = 0.0; if(fator > 0.5) fator = 0.5; if(tipo == BOLA) { if(canal == RED) valor = MULT*(fator)*cena->bolas[indice].superficie.difusa.c[0]; if(canal == GREEN) valor = MULT*(fator)*cena->bolas[indice].superficie.difusa.c[1]; if(canal == BLUE) valor = MULT*(fator)*cena->bolas[indice].superficie.difusa.c[2]; } else { if(canal == RED) valor = MULT*(fator)*cubos[indice].superficie.difusa.c[0]; if(canal == GREEN) valor = MULT*(fator)*cubos[indice].superficie.difusa.c[1]; if(canal == BLUE) valor = MULT*(fator)*cubos[indice].superficie.difusa.c[2]; } } else { valor = 0.0; } } else { valor = 0.0; } return valor; } else { return 0.0; } } //radiosidade.h MatrizEsparsa *ObtemMatrizDeTransferenciaDasLuzes(Cena *cena, int canal, BaseAprox *bas); /* Obtém a matriz de transferencia de radiância {D} entre as luzes externas e os sítios de amostragem {bas->sd[0..bas->nel-1]}. Se {bas->Dd} for {NULL}, constrói a matriz com {CalculaMatrizDeTransferenciaDasLuzes}, e salva o resultado em {bas->Dd}. Em qualquer caso, devolve {bas->Dd}. */ double ShepardPesos(Ponto3d P, int N, Cena *cena, Vetor3d Normal,int canal) { int i = 0; double somatorio = 0.0; //somatorio dos pesos Vetor3d ni; double ang = 0.0; double d_p = 0.0; for(i = 0; i < N; i++) { //encontra a normal do s i if(bas->sb[i].tipo == CUBO) EncontraNormalCubo1(cubos[bas,Vis,Cam.ind[i]],bas->sb[i],&ni,bas,Vis,Cam.l[i]); if(bas->sb[i].tipo == BOLA) EncontraNormalBola(cena->bolas[bas,Vis,Cam.ind[i]],bas->sb[i],&ni); //normaliza as normais NormalizaVetor(&ni); NormalizaVetor(&Normal); //calcula o cos do angulo ang = (ProdutoEscalar(ni,Normal)); if(ang <= 0.00001) ang = DELTA; //base de Shepard if(Base == SHEPARD) CoeficienteAux[i]= (1/(pow((Distancia2Pontos(P,bas->sb[i]))/ang,4))); else //base do artigo meshless if(Base == RADIAL_SHEPARD) { double dist = Distancia2Pontos(P,bas->sb[i])+ALPHA; double d = ((1.0)/sqrt(Alpha[i])); d_p = dist/d; if(d_p <= 1.0) dist = 2*pow(d_p,3) - 3*pow(d_p,2) + 1; else dist = ALPHA; CoeficienteAux[i] = dist*ang; //printf(""); } somatorio += CoeficienteAux[i]; } for(i = 0; i < N; i++) { CoeficienteAux[i] = (CoeficienteAux[i]/somatorio); } return somatorio; } void InterpolaRadianciaDaBaseShepard(Ponto3d P, int N, frgb_t *cor,Cena *cena, Vetor3d Normal) { double somatorio = 0.0; int i = 0; somatorio = ShepardPesos(P,N,bolas,cubos,Normal,TODAS); cor->r = 0.0; cor->g = 0.0; cor->b = 0.0; for(i = 0; i < N; i++) { //if(i == INDICE1 || i == INDICE2 || i == INDICE3) { cor->r += (CoeficienteAux[i])*bas,Vis,Cam.c[i].c[0]; cor->g += (CoeficienteAux[i])*bas,Vis,Cam.c[i].c[1]; cor->b += (CoeficienteAux[i])*bas,Vis,Cam.c[i].c[2]; } } } void ReconstrucaoShepard(Filme *filme, int N, Cena *cena, int num_Pontos) { int i = 0; int indice = 0; frgb_t c; Vetor3d ni; for(i = 0; i < N ; i++) { if(bas,Vis,Cam.tipo_i[i] == CUBO) EncontraNormalCubo1(cubos[bas,Vis,Cam.ind_i[i]],bas,Vis,Cam.s_i[i],&ni,bas,Vis,Cam.l_i[i]); if(bas,Vis,Cam.tipo_i[i] == BOLA) EncontraNormalBola(cena->bolas[bas,Vis,Cam.ind_i[i]],bas,Vis,Cam.s_i[i],&ni); //normaliza a normal NormalizaVetor(&ni); InterpolaRadianciaDaBaseShepard(bas,Vis,Cam.s_i[i],bas->nel,&c,bolas,cubos,ni);//reconstroi o ponto i indice = bas,Vis,Cam.h_pixel[i] + (filme->n_colunas)*bas,Vis,Cam.v_pixel[i];//calcula o indice do pixel filme->pixel[indice].cor = c;//"pinta" o pixel da filme } for(i = 0; i < bas->nel ; i++) { indice = bas,Vis,Cam.x[i] + (filme->n_colunas)*bas,Vis,Cam.y[i];//calcula o indice do pixel filme->pixel[indice].cor = bas,Vis,Cam.c[i];//"pinta" o pixel da filme } } void CalculaMatrizesSis,Vis,Cam(int N, int canal, Cena *cena) { int i,j; Vetor3d nj; double *a; a = (double*)malloc(sizeof(double)*N*N); if(a == NULL) { printf("\n ERRO DE ALOCACAO - Funcao:CalculaMatrizesSis,Vis,Cam"); exit(1); } for(j = 0; j < N; j++) { //encontra normal do centro j if(bas->el[j].centro.tipo == CUBO) EncontraNormalCubo1(cubos[bas,Vis,Cam.ind[j]],bas->el[j].centro,&nj,bas,Vis,Cam.l[j]); if(bas->el[j].centro.tipo == BOLA) EncontraNormalBola(bolas[bas,Vis,Cam.ind[j]],bas->el[j].centro,&nj); if(Base == RADIAL_SHEPARD) { ShepardPesos(bas->el[j].centro,bas->nel,bolas,cubos,nj,TODAS); for(i = 0; i < N; i++) { a[j*N+i] = CoeficienteAux[i]; } } if(Base == RADIAL) { for(i = 0; i < N; i++) { Vetor3d ni; if(bas->sb[i].tipo == CUBO) EncontraNormalCubo1(cubos[bas,Vis,Cam.ind[i]],bas->sb[i],&ni,bas,Vis,Cam.l[i]); if(bas->sb[i].tipo == BOLA) EncontraNormalBola(cena->bolas[bas,Vis,Cam.ind[i]],bas->sb[i],&ni); a[j*N+i] = RBF_Gaussiana(bas->el[j].centro,bas->sb[i],nj,ni,j); } } } MatrizParaMatrizEsparsaAloca(a,N,N,&MM); free(a); a = NULL; } double RBF_Gaussiana(Ponto3d A, Ponto3d B, Vetor3d Na, Vetor3d Nb, int Ind) { double dist = Distancia2Pontos(A,B); double d = ((1.0)/sqrt(Alpha[Ind])); double ang = 0.0; //normaliza vetores NormalizaVetor(&Na); NormalizaVetor(&Nb); ang = (ProdutoEscalar(Na,Nb)); if(ang <= 0.01) ang = 0.000000000000000000000000001; /*double d_p = 0.0; d_p = dist/d; if(d_p <= 1.0) dist = 2*pow(d_p,3) - 3*pow(d_p,2) + 1; else dist = 0.0; dist = dist*ang; return dist; */ if(dist > 2.0*d)// || ang < 0.5) return(0.0); dist = dist/ang; //Gausianas return(exp(-3.0*Alpha[Ind]*(dist*dist)) ); //return(exp(-1.0*Alpha[Ind]*(dist*dist*(d/1.5)))); } //geometria_da_cena.h //dados de uma fonte de luz exerna: typedef struct Luz { Ponto3d posicao; //posicao da fonte no espaço. double brilho; //intensidade nominal da fonte. double d_ref; //distancia onde a radiância máxima é {brilho}. } Luz; // !!! as luzes são desnecessárias -- basta usar objetos luminosos. int num_luzes; //número de luzes externas. Luz *luzes; //as luzes, indexadas de 0 a {num_luzes-1}. //radiosidade.c void AlocaBaseAprox(BaseAprox *bas, int max_nel); /* Aloca vetores internos de {*bas} com espaço para {max_nel} elementos. */ void ExpandeBaseAprox(BaseAprox *bas, int indice); /* Realoca a base {bas} de modo que exista o elemento {bas->el[indice]}. */ void AlocaBaseAprox(BaseAprox *bas, int max_nel) { //!!! deveria alocar também o cabeçalho. bas->max_nel = max_nel; bas->el = (ElemFinito *)NOTNULL(malloc(sizeof(ElemFinito)*max_nel)); fprintf(stderr, "alocados %d elementos para a base\n", max_nel); } void ExpandeBaseAprox(BaseAprox *bas, int indice) { if (indice >= bas->max_nel) { //calcula novo tamanho {novo_max}, {50%} maior que o velho: int novo_max = 3*bas->max_nel/2; //garante que é suficiente para o elemento {indice}: if (indice >= novo_max) { novo_max = indice + 1; } novo_max += 1000; //para garantir um aumento mínimo. //realoca vetor e lembra o novo tamanho: bas->el = (ElemFinito *)NOTNULL(realloc(bas->el, novo_max*sizeof(ElemFinito))); bas->max_nel = novo_max; } } //matriz-esparsa.c void InverteMatrizEsparsa(MatrizEsparsa *A, MatrizEsparsa *B); /** InverteMatrizEsparsa: Funcao que recebe uma matriz esparsa, realiza a inversao -------------------- e joga o resultado em outra matriz esparsa **/ void InverteMatrizEsparsa(Matriz_Esparsa *A, Matriz_Esparsa *B) { int i = 0; double *a; double *inv; double *res; a = (double*)malloc(sizeof(double)*A->n_linhas*A->n_colunas); inv = (double*)malloc(sizeof(double)*A->n_linhas*A->n_colunas); res = (double*)malloc(sizeof(double)*A->n_linhas*A->n_colunas); if(a == NULL || inv == NULL || res == NULL) { printf("\n ERRO DE ALOCACAO - Funcao:InverteMatrizEsparsa"); exit(1); } for(i = 0; i < A->n_linhas*A->n_colunas; i++) inv[i] = a[i] = 0.0; for(i = 0; i < A->num_elementos; i++) { a[A->a[i].i*Sistema.n_centros + A->a[i].j] = A->a[i].aij; } invert_matrix (a,inv,Sistema.n_centros); MultiplicaMatrizMatriz(a,inv,res,Sistema.n_centros); MatrizParaMatrizEsparsaAloca(inv,Sistema.n_centros,Sistema.n_centros,B); free(a); free(inv); free(res); a = inv = res = NULL; } /** InsereColunaMatrizEsparsa:Funcao que insere uma Coluna de uma matriz esparsa ------------------------- **/ void InsereColunaMatrizEsparsa(Matriz_Esparsa *M, double *C, int Coluna) { int i; for(i = 0; i < M->n_linhas; i++) { if(C[i] != 0.0) { InsereElementoMatrizEsparsa(M,C[i],i,Coluna); } } } /** RetornaColunaMatrizEsparsa:Funcao que retorna uma coluna de uma matriz esparsa ------------------------- **/ void RetornaColunaMatrizEsparsa(Matriz_Esparsa *M, double *C, int Coluna) { int i; for(i = 0; i < M->n_linhas; i++) C[i] = 0.0; for(i = 0; i < M->num_elementos; i++) { //se for a coluna desejada if(M->a[i].j == Coluna) { C[M->a[i].i] = M->a[i].aij; } } } /** RetornaLinhaMatrizEsparsa:Funcao que retorna uma linha de uma matriz esparsa ------------------------- **/ void RetornaLinhaMatrizEsparsa(Matriz_Esparsa *M, double *L, int Linha) { int i; for(i = 0; i < M->num_colunas; i++) L[i] = 0.0; for(i = 0; i < M->num_elementos; i++) { //se for alinha desejada if(M->a[i].i == Linha) { L[M->a[i].j] = M->a[i].aij; } //if(M->a[i].i > Linha) //break; } } /** MatrizParaMatrizEsparsa: Funcao que "converte" uma matriz convencional em ----------------------- uma matriz esparsa **/ /** ContaValoresNaoNulos: Funcao que conta o numero de elementos nao -------------------- nulos de uma matriz **/ int ContaValoresNaoNulos(double *M, int N_Linhas, int Num_Colunas) { int cont = 0; int i,j; for(i = 0; i < N_Linhas; i++) { for(j = 0; j < Num_Colunas; j++) { if(M[i*Num_Colunas + j] != 0.0) cont++; } } return cont; } /** InsereElementoMatrizEsparsa : Funcao que insere um elemento em uma --------------------------- matriz esparsa **/ void InsereElementoMatrizEsparsa(Matriz_Esparsa *M, double Val, int I, int J) { M->a[M->num_elementos].i = I; M->a[M->num_elementos].j = J; M->a[M->num_elementos].aij = Val; M->num_elementos++; } /** DeslocaMatrizEsparsa: Funcao que desaloca uma matriz esparsa ------------------- **/ void DesalocaMatrizEsparsa(Matriz_Esparsa *M) { free(M->a); M->a = NULL; printf("\n Desalocando.... (%dx%d) e %d elementos ",M->n_linhas,M->num_colunas,M->num_elementos); M->n_linhas = M->num_colunas = M->num_elementos = 0; } /** AlocaMatrizEsparsa: Funcao que aloca uma matriz esparsa ------------------- **/ void AlocaMatrizEsparsa(Matriz_Esparsa *M, int Num, int N_Linhas,int Num_Colunas) { M->a = (Elemento*)malloc(sizeof(Elemento)*(Num)); M->num_elementos = 0; M->n_linhas = N_Linhas; M->num_colunas = Num_Colunas; return; } //geometria_da_cena.c void DeterminaNormalDeCubo(Cubo *obj, Ponto3d *posicao, Vetor3d *normal) { //inicializa o vetor normal normal->c[0] = 0.0; normal->c[1] = 0.0; normal->c[2] = 0.0; switch(Li) { case(0): normal->c[0] = 1.0; break; case(1): normal->c[0] = -1.0; break; case(2): normal->c[1] = 1.0; break; case(3): normal->c[1] = -1.0; break; case(4): normal->c[2] = 1.0; break; case(5): normal->c[2] = -1.0; break; } } /** EncontraNormalCubo: Funcao que dado um Cubo e o ponto de intersecção retorna o vetor normal **/ void EncontraNormalCubo(Cubo cubo, Ponto3d sint, Vetor3d *normal) { //inicializa o vetor normal normal->c[0] = 0.0; normal->c[1] = 0.0; normal->c[2] = 0.0; switch(cubo.li) { case(0): normal->c[0] = 1.0; break; case(1): normal->c[0] = -1.0; break; case(2): normal->c[1] = 1.0; break; case(3): normal->c[1] = -1.0; break; case(4): normal->c[2] = 1.0; break; case(5): normal->c[2] = -1.0; break; } } /** EncontraNormalInterseccao: Funcao que dado o ponto de intersecção e informacoes da esfera retorna o vetor normal a superficie **/ void EncontraNormalInterseccao(Bola obj, Ponto3d cam->olho, Vetor3d dir, double t, Vetor3d *normal,Ponto3d *sint) { double result; //encontra o ponto de sint r3_dir(&dir); sint->c[0] = cam->olho.c[0] + t*dir.c[0]; sint->c[1] = cam->olho.c[1] + t*dir.c[1]; sint->c[2] = cam->olho.c[2] + t*dir.c[2]; //encontra a normal a superficie normal->c[0] = (sint->c[0] - obj.c[0])/obj.raio; normal->c[1] = (sint->c[1] - obj.c[1])/obj.raio; normal->c[2] = (sint->c[2] - obj.c[2])/obj.raio; result = pow((sint->c[0] -obj.c[0]),2) + pow((sint->c[1] -obj.c[1]),2) + pow((sint->c[2] -obj.c[2]),2) - pow(obj.raio,2); } /** RetornaPontonInterseccao: Funcao que dado os parametros retorna o ponto de intersecção **/ Ponto3d RetornaPontonInterseccao(Ponto3d cam->olho, Vetor3d dir, double t) { Ponto3d sint; r3_dir(&dir); sint.c[0] = cam->olho.c[0] + t*dir.c[0]; sint.c[1] = cam->olho.c[1] + t*dir.c[1]; sint.c[2] = cam->olho.c[2] + t*dir.c[2]; return(sint); } MatrizEsparsa *ObtemMatrizDeTransferenciaDosCentros(Opcoes *o, Cena *cena, int C, VetorDeSitios *st); /* Estas funções obtém a matriz {T} de transferência de radiância entre os sítios {st} para o canal {C} e para a {cena} dada. Se {PR_ENT = o->prefixoDasMatsDeEntrada} não é {NULL}, lê a matriz do arquivo "{PR_ENT}_T.mat", senão calcula a matriz. Em qualquer caso, se {PR_SAI = o->prefixoSaida} não é NULL, grava a matriz no arquivo "{PR_SAI}_T.mat". */ //geometria_basica.c void InverteVetor(Vetor3d *v); double Distancia2Pontos(Ponto3d P1, Ponto3d P2); void ProdutoVetorial(Vetor3d A, Vetor3d B, Vetor3d *C); double ProdutoEscalar(Vetor3d A, Vetor3d B); void InverteVetor(Vetor3d *v) { v->c[0] = -v->c[0]; v->c[1] = -v->c[1]; v->c[3] = -v->c[3]; } double Distancia2Pontos(Ponto3d P1, Ponto3d P2) { double dist = 0.0; //calcula a distancia entre os pontos dist = pow(P1.c[0] - P2.c[0],2.0) + pow(P1.c[1] - P2.c[1],2.0) + pow(P1.c[3] - P2.c[3],2.0); dist = sqrt(dist); return (dist); } void ProdutoVetorial(Vetor3d A, Vetor3d B, Vetor3d *C) { C->c[0] = A.c[1]*B.c[3] - A.c[3]*B.c[1]; C->c[1] = A.c[3]*B.c[0] -A.c[0]*B.c[3]; C->c[3] = A.c[0]*B.c[1] - A.c[1]*B.c[0]; return; } double ProdutoEscalar(Vetor3d A, Vetor3d B) { double produto_escalar = 0.0; produto_escalar = A.c[0]*B.c[0] + A.c[1]*B.c[1] + A.c[3]*B.c[3]; return(produto_escalar); } Vetor3d DirecaoAleatoria(void) { //!!! nao é uniforme na esfera. Vetor3d v; v.c[0] = random(); v.c[1] = random(); v.c[3] = random(); r3_dir_safe(&v); return v; } //radioso.c //sítios visíveis na cena correspondentes aos pixeis da imagem. typedef struct SitiosDePixeis { int max_nvs; //tamanho dos vetores. int nvs; //número de sítios no conjunto. Sitio *sitio_i; //sitios dos pixeis, indexados 0 a {nvs-1}. int *h_pixel; //coordenada do pixel int *v_pixel; //coordenada do pixel } SitiosDePixeis; void AlocaSitiosDePixeis(SitiosDePixeis *vis, int max_nvs); /* Aloca vetores internos de {*vis} com espaço para {max_nvs} sítios. */ void ExpandeSitiosDePixeis(SitiosDePixeis *vis, int indice); /* Realoca os vetores internos de {vis} de modo que exista o sítio {vis->sitio_i[indice]}. */ void DesalocaSitiosDePixeis(SitiosDePixeis *vis); /* Libera a memória usada pelo conjunto de sítios {vis} (menos o cabeçalho). */ printf("\n Tempo de execucao: \n"); printf(" Execucao Rendering: %lf s \n",difftime(rend_fim,inicio)); printf(" Execucao Precomputacao: %lf s \n",difftime(pre_fim,pre_ini)); printf(" Execucao Multiplicacao: %lf s \n",difftime(mul_fim,mul_ini)); printf(" Execucao sis,Vis,Cam: %lf s \n",difftime(sis_fim,sis_ini)); printf(" Execucao Imagem: %lf s \n",difftime(img_fim,img_ini)); printf(" EXECUCAO TOTAL: %lf s \n",difftime(fim,inicio)); /** AlocaPrecomputacao: Funcao que realiza a alocação das variaveis de precomputacao **/ void AlocaPrecomputacao(BaseAprox *sis) { Ir = (double*)malloc(sizeof(double)*(sis->nel)); Ig = (double*)malloc(sizeof(double)*(sis->nel)); Ib = (double*)malloc(sizeof(double)*(sis->nel)); if(Ir == NULL || Ig == NULL || Ib == NULL) { printf("\n ERRO ALOCACAO - Funcao:AlocaPrecomputacao"); exit(1); } } /** DesalocaPrecomputacao: Funcao que realiza a desalocação das variaveis de precomputacao **/ void DesalocaPrecomputacao(void) { free(Ir) ; Ir = NULL; free(Ig); Ig = NULL; free(Ib); Ib = NULL; } /** AlocaVariaveis: Funcao que realiza a alocação das variaveis **/ void AlocaVariaveis(BaseAprox *sis) { CoeficienteRed = (double*)NOTNULL(malloc(sizeof(double)*(sis->nel))); CoeficienteGreen = (double*)NOTNULL(malloc(sizeof(double)*(sis->nel))); CoeficienteBlue = (double*)NOTNULL(malloc(sizeof(double)*(sis->nel))); CoeficienteAux = (double*)NOTNULL(malloc(sizeof(double)*(sis->nel))); } /** DesalocaVariaveis: Funcao que realiza a desalocação das variaveis **/ void DesalocaVariaveis(void) { free(CoeficienteRed); CoeficienteRed = NULL; free(CoeficienteGreen); CoeficienteGreen = NULL; free(CoeficienteBlue); CoeficienteBlue = NULL; free(CoeficienteAux); CoeficienteAux = NULL; } void ExpandeSitiosDePixeis(SitiosDePixeis *vis, int indice) { if (indice >= vis->max_nvs) { //calcula novo tamanho {novo_max}, {50%} maior que o velho: int novo_max = 3*vis->max_nvs/2; //garante que é suficiente para o elemento {indice}: if (indice >= novo_max) { novo_max = indice + 1; } novo_max += 1000; //para garantir um aumento mínimo. //realoca vetor e lembra o novo tamanho: vis->sitio_i = (Sitio *)NOTNULL(realloc(vis->sitio_i, novo_max*sizeof(Sitio))); vis->h_pixel = (int *)NOTNULL(realloc(vis->h_pixel, novo_max*sizeof(int))); vis->v_pixel = (int *)NOTNULL(realloc(vis->v_pixel, novo_max*sizeof(int))); vis->max_nvs = novo_max; } } void DesalocaSitiosDePixeis(SitiosDePixeis *vis) { free(vis->sitio_i); vis->sitio_i = NULL; free(vis->h_pixel); vis->h_pixel = NULL; free(vis->v_pixel); vis->v_pixel = NULL; } /* ---------------------------------------------------------------------- */ //encontra a normal do centro i if(sis.el[i].centro.tipo == CUBO) EncontraNormalCubo1(cubos[sis,Vis,Cam.ind[i]],sis.el[i].centro,&ni,sis,Vis,Cam.l[i]); if(sis.el[i].centro.tipo == BOLA) EncontraNormalBola(cena->bolas[sis,Vis,Cam.ind[i]],sis.el[i].centro,&ni); NormalizaVetor(&ni); /* ---------------------------------------------------------------------- */ //encontra a normal do centro j if(sis.el[j].centro.tipo == CUBO) EncontraNormalCubo1(cubos[sis,Vis,Cam.ind[j]],sis.el[j].centro,&nj,sis,Vis,Cam.l[j]); if(sis.el[j].centro.tipo == BOLA) EncontraNormalBola(cena->bolas[sis,Vis,Cam.ind[j]],sis.el[j].centro,&nj); NormalizaVetor(&nj); //distancia entre as normais double ang = (ProdutoEscalar(ni,nj)); if(ang < DELTA) { ang = DELTA; } double dist = d/ang; // Distancia dos sítios. /* ---------------------------------------------------------------------- */ if(tipo == BOLA) { if(C == RED) Mat[j*num_luzes + i] = (fator)*cena->bolas[indice].superficie.difusa.c[0]; if(C == GREEN) Mat[j*num_luzes + i] = (fator)*cena->bolas[indice].superficie.difusa.c[1]; if(C == BLUE) Mat[j*num_luzes + i] = (fator)*cena->bolas[indice].superficie.difusa.c[2]; } else { if(C == RED) Mat[j*num_luzes + i] = (fator)*cubos[indice].superficie.difusa.c[0]; if(C == GREEN) Mat[j*num_luzes + i] = (fator)*cubos[indice].superficie.difusa.c[1]; if(C == BLUE) Mat[j*num_luzes + i] = (fator)*cubos[indice].superficie.difusa.c[2]; } } else { Mat[j*num_luzes + i] = 0.0; } } else { Mat[j*num_luzes + i] = 0.0; } /* ---------------------------------------------------------------------- */ void InsereDeslocaVetor(double *V,int Pos,double Val); /* Insere um valor na posicao Pos e desloca os valores restantes em 1 unidade. */ void InsereDeslocaVetor(double *V,int Pos,double Val) { int i = 0; double origem, destino; origem = V[Pos]; for(i = Pos; i < LIM; i++) { destino = V[i+1]; V[i+1] = origem; origem = destino; } V[Pos] = Val; }