#include #include #include #include #include /*Validation functions*/ int NoValidation(double* c, int basis_size, bool_t* validity_array){ return basis_size; } /*Argpareser reader*/ lighting_data_model_t parse_lighting_model(argparser_t* pp, bool_t determineParams){ lighting_data_model_t lmd; lmd.pointlikeModel = POINTLIKE_MODEL_PLAIN; if (argparser_keyword_present_next(pp, "pointlike")){ lmd.lightingModelType = POINTLIKE_MODEL; fprintf(stderr," POINTLIKE_MODEL\n"); if(determineParams){ if(argparser_keyword_present_next(pp, "backplane")){ lmd.pointlikeModel = POINTLIKE_MODEL_BACKP; fprintf(stderr," Using backplane : %d\n",lmd.pointlikeModel); } if(argparser_keyword_present_next(pp, "pure")){ lmd.pointlikeModel = POINTLIKE_MODEL_LPURE; fprintf(stderr," Using pure model : %d\n",lmd.pointlikeModel); } lmd.pointlikeK = INF; lmd.pointlikeRho = INF; if(argparser_keyword_present_next(pp, "highlightK")){ lmd.pointlikeK = argparser_get_next_double(pp, 0.0, INF); fprintf(stderr," Using highlight sharpness: %lf\n",lmd.pointlikeK); } if(argparser_keyword_present_next(pp, "terminatorRho")){ lmd.pointlikeRho = argparser_get_next_double(pp, 0.0, INF); fprintf(stderr," Using terminator sharpness: %lf\n",lmd.pointlikeRho); } } }else if (argparser_keyword_present_next(pp, "harmonic")){ lmd.lightingModelType = HARMONIC_MODEL; fprintf(stderr," HARMONIC_MODEL\n"); if(determineParams){ argparser_get_keyword_next(pp, "degree"); lmd.harmonicDegree = argparser_get_next_int(pp, 1, 10000); fprintf(stderr," degree %d\n",lmd.harmonicDegree); } }else if (argparser_keyword_present_next(pp, "radial")){ lmd.lightingModelType = RADIALBASIS_MODEL; fprintf(stderr," RADIALBASIS_MODEL\n"); if(determineParams){ argparser_get_keyword_next(pp, "resolution"); lmd.radialResolution = argparser_get_next_int(pp, 1, 10000); argparser_get_keyword_next(pp, "span"); lmd.radialSpan =argparser_get_next_double(pp, 0.0, 10000.0); fprintf(stderr," resolution %d\n span %4.4lf\n",lmd.radialResolution,lmd.radialSpan); } }else if (argparser_keyword_present_next(pp, "glossy")){ lmd.lightingModelType = GLOSSYLIKE_MODEL; fprintf(stderr," GLOSSYLIKE_MODEL\n"); if(determineParams){ argparser_get_keyword_next(pp, "glossiness"); lmd.glossylikeG0 = argparser_get_next_double(pp, 0.0, 10000); lmd.glossylikeG1 = argparser_get_next_double(pp, 0.0, 10000); fprintf(stderr," G0 %4.4lf\n G1 %4.4lf\n",lmd.glossylikeG0,lmd.glossylikeG1); } }else if (argparser_keyword_present_next(pp, "harmonicSG")){ lmd.lightingModelType = HARMONICSG_MODEL; fprintf(stderr," HARMONICSG_MODEL\n"); if(determineParams){ argparser_get_keyword_next(pp, "degree"); lmd.harmonicDegree = argparser_get_next_int(pp, 1, 10000); fprintf(stderr," degree %d\n",lmd.harmonicDegree); } }else if (argparser_keyword_present_next(pp, "compact")){ lmd.lightingModelType = COMPACT_MODEL; fprintf(stderr," COMPACT_MODEL\n"); if(determineParams){ argparser_get_keyword_next(pp, "dirlight"); lmd.compactU.c[0] = argparser_get_next_double(pp, 0.0, 10000); lmd.compactU.c[1] = argparser_get_next_double(pp, 0.0, 10000); lmd.compactU.c[2] = argparser_get_next_double(pp, 0.0, 10000); r3_t normU; double lenghtU = r3_dir(&(lmd.compactU),&normU); if(lenghtU != 1.0){ fprintf(stderr,"Non normalized U (%9.6lf %9.6lf %9.6lf), normalized to (%9.6lf %9.6lf %9.6lf)\n",lmd.compactU.c[0],lmd.compactU.c[1],lmd.compactU.c[2],normU.c[0],normU.c[1],normU.c[2]); } lmd.compactU = normU; argparser_get_keyword_next(pp, "dirbackplane"); lmd.compactXi.c[0] = argparser_get_next_double(pp, 0.0, 10000); lmd.compactXi.c[1] = argparser_get_next_double(pp, 0.0, 10000); lmd.compactXi.c[2] = argparser_get_next_double(pp, 0.0, 10000); r3_t normXi; double lenghtXi = r3_dir(&(lmd.compactXi),&normXi); if(lenghtXi != 1.0){ fprintf(stderr,"Non normalized Xi (%9.6lf %9.6lf %9.6lf), normalized to (%9.6lf %9.6lf %9.6lf)\n",lmd.compactXi.c[0],lmd.compactXi.c[1],lmd.compactXi.c[2],normXi.c[0],normXi.c[1],normXi.c[2]); } argparser_get_keyword_next(pp, "rho"); lmd.compactRho = argparser_get_next_double(pp, 0.0, 10000); argparser_get_keyword_next(pp, "K"); lmd.compactK = argparser_get_next_double(pp, 0.0, 10000); } }else{ fprintf(stderr,"Keyword was %s\n",argparser_get_next(pp)); demand(FALSE, "unknown lighting model"); } return lmd; } /*Generic model type*/ approx_model_t* create_approx_lighting_model(model_type_t type){ approx_model_t* am = (approx_model_t*) malloc(sizeof(approx_model_t)); if(type == HARMONIC_MODEL){ am->phi = harmonic_phi; am->get_num_components = harmonic_get_number_components; am->set_alphas = harmonic_set_alphas; am->get_alphas = harmonic_get_alphas; am->validade_results = NULL; am->copy_data = harmonic_copy_lighting_data; am->release_data = harmonic_release_lighting_data; am->evaluate = harmonic_shading; am->write_params = harmonic_write_params; am->read_params = harmonic_read_params; /*Non linear */ am->update_weights = NULL; am->compare = NULL; am->get_num_nl_params = NULL; am->pack_nl_params = NULL; am->unpack_nl_params = NULL; am->pack_nl_param_errors = NULL; am->unpack_nl_param_errors = NULL; }else if(type == RADIALBASIS_MODEL){ am->phi = radial_phi; am->get_num_components = radial_get_number_components; am->set_alphas = radial_set_alphas; am->get_alphas = radial_get_alphas; am->validade_results = NULL; am->copy_data = radial_copy_lighting_data; am->release_data = radial_release_lighting_data; am->evaluate = radial_shading; am->write_params = radial_write_params; am->read_params = radial_read_params; /*Non linear */ am->update_weights = NULL; am->compare = NULL; am->get_num_nl_params = NULL; am->pack_nl_params = NULL; am->unpack_nl_params = NULL; am->pack_nl_param_errors = NULL; am->unpack_nl_param_errors = NULL; }else if(type == HARMONICSG_MODEL){ am->phi = harmonicSG_phi; am->get_num_components = harmonicSG_get_number_components; am->set_alphas = harmonicSG_set_alphas; am->get_alphas = harmonicSG_get_alphas; am->validade_results = NULL; am->copy_data = harmonicSG_copy_lighting_data; am->release_data = harmonicSG_release_lighting_data; am->evaluate = harmonicSG_shading; am->write_params = harmonicSG_write_params; am->read_params = harmonicSG_read_params; /*Non linear */ am->update_weights = NULL; am->compare = NULL; am->get_num_nl_params = NULL; am->pack_nl_params = NULL; am->unpack_nl_params = NULL; am->pack_nl_param_errors = NULL; am->unpack_nl_param_errors = NULL; }else if( type == COMPACT_MODEL){ am->phi = compact_phi; am->get_num_components = compact_get_number_components; am->set_alphas = compact_set_alphas; am->get_alphas = compact_get_alphas; am->validade_results = NULL; am->copy_data = compact_copy_lighting_data; am->release_data = compact_release_lighting_data; am->evaluate = compact_shading; am->write_params = compact_write_params; am->read_params = compact_read_params; /*Non linear */ am->update_weights = compact_update_weights; am->compare = compact_compare_lightings; am->get_num_nl_params = compact_num_nl_params; am->pack_nl_params = compact_pack_parameter; am->unpack_nl_params = compact_unpack_parameter; am->pack_nl_param_errors = NULL; am->unpack_nl_param_errors = NULL; }else{ free(lm); return NULL; } return lm; } /*Polynomial model*/ harmonic_lighting_t* harmonic_init_components(int g,r3_t view_dir){ harmonic_lighting_t* pl = (harmonic_lighting_t*)malloc(sizeof(harmonic_lighting_t)); int num_comp = (g+1)*(g+1); pl->num_comp = num_comp; pl->view_dir = view_dir; pl->coeficients = (double*)malloc(sizeof(double)*num_comp); pl->triplets = (double**)malloc(sizeof(double*)*num_comp); int i,j,k; for(i = 0; i < num_comp; i++){ pl->coeficients[i] = 1.0; pl->triplets[i] = (double*)malloc(sizeof(double)*3); } int count = 0; for(i = g; i >= 0; i--){ for( j = g-i; j >= 0; j--){ k = g - i - j; pl->triplets[count][0] = i; pl->triplets[count][1] = j; pl->triplets[count][2] = k; count++; } } g = g -1; for(i = g; i >= 0; i--){ for( j = g-i; j >= 0; j--){ k = g - i - j; pl->triplets[count][0] = i; pl->triplets[count][1] = j; pl->triplets[count][2] = k; count++; } } pl->num_comp = count; return pl; } double harmonic_phi(int r, double* x, void* l_data){ r3_t normal = (r3_t){{x[0],x[1],x[2]}}; harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; int ri = pl->triplets[r][0]; int rj = pl->triplets[r][1]; int rk = pl->triplets[r][2]; double phiR = (pow(normal.c[0],ri))*(pow(normal.c[1],rj))*(pow(normal.c[2],rk)); //phiR = phiR*pl->weights[r]; return phiR; } int harmonic_get_number_components(void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; return pl->num_comp; } void harmonic_set_alphas(double* C,void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; int i; for(i = 0; i < n; i++){ pl->coeficients[i] = C[i]; } } void harmonic_get_alphas(void* l_data,double* C,int n){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; assert(n == pl->num_comp); int i; for(i = 0; i < pl->num_comp; i++){ C[i] = pl->coeficients[i]; } } void* harmonic_copy_lighting_data(void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; harmonic_lighting_t* npl = (harmonic_lighting_t*) malloc(sizeof(harmonic_lighting_t)); npl->num_comp = pl->num_comp; npl->triplets = (double**)malloc(sizeof(double*)*(npl->num_comp)); npl->coeficients = (double*)malloc(sizeof(double)*(npl->num_comp)); npl->view_dir = pl->view_dir; int i; for(i = 0; i < npl->num_comp;i++){ npl->coeficients[i] = pl->coeficients[i]; npl->triplets[i] = (double*)malloc(sizeof(double)*3); int j; for(j = 0; j < 3;j++){ npl->triplets[i][j] = pl->triplets[i][j]; } } return npl; } void harmonic_release_lighting_data(void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; int i; for(i = 0; i < pl->num_comp;i++){ free(pl->triplets[i]); } free(pl->triplets); free(pl->coeficients); free(l_data); } double harmonic_shading(double* x,void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; int r; double val = 0; for(r = 0; r < pl->num_comp; r++){ val+= (pl->coeficients[r]*harmonic_phi(r,x,l_data)); } return val; } void harmonic_write_params(FILE* arq,void* l_data){ harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data; fprintf(arq," %d\n",pl->num_comp); int i; for(i = 0; i < pl->num_comp;i++){ fprintf(arq, "%8.6lf %8.6lf %8.6lf %8.6lf\n",pl->coeficients[i],pl->triplets[i][0],pl->triplets[i][1],pl->triplets[i][2] ); } } void* harmonic_read_params(FILE* arq){ harmonic_lighting_t* pl = (harmonic_lighting_t*) malloc(sizeof(harmonic_lighting_t)); fscanf(arq,"%d",&(pl->num_comp)); int i; int num_comp = pl->num_comp; pl->coeficients = (double*)malloc(sizeof(double)*num_comp); pl->triplets = (double**)malloc(sizeof(double*)*num_comp); for(i = 0; i < pl->num_comp;i++){ pl->triplets[i] = (double*)malloc(sizeof(double)*3); fscanf(arq, "%lf %lf %lf %lf",&(pl->coeficients[i]),&(pl->triplets[i][0]),&(pl->triplets[i][1]),&(pl->triplets[i][2]) ); } return pl; } void harmonic_update_weights(void* l_data,double** X, double* F,double* weights, int n,double * sigma){ double sumd2 = 0.0; double sumw = 1.0e-200; double sigma2 = *sigma; sigma2 = sigma2*sigma2; int i; for(i = 0; i< n; i++){ // r3_t normal = (r3_t){{X[i][0],X[i][1],X[i][2]}}; double computedGO = harmonic_shading(X[i],l_data); double originalGO = F[i]; double d2 = (computedGO-originalGO)*(computedGO-originalGO); weights[i] = (d2 > 5.0*sigma2 ? 0.0 : exp(-d2/(2.0*sigma2)) ); sumd2+= d2; sumw+=weights[i]; } *sigma = sqrt(sumd2/(double)n); } double harmonic_compare_lightings(void* l_data1,void* l_data2){ return 0.0; } /*Radial Basis Model*/ radialbasis_lighting_t* radialbasis_init_components(int resolution,double span,r2_t center,double raio, r2_t estica,r3_t pole){ r2_t* pontos; r3_t* normals; int num_comp; gera_pontos_no_gabarito_eliptico(resolution,center,raio,estica,&pontos, &normals,&num_comp,pole); // gera_normais_no_gabarito(resolution,center, num_comp,pole); radialbasis_lighting_t* rl = (radialbasis_lighting_t*)malloc(sizeof(radialbasis_lighting_t)); rl->num_comp = num_comp; rl->bradius = (double*)malloc(sizeof(double)*num_comp); rl->coeficients = (double*)malloc(sizeof(double)*num_comp); rl->bcenter = normals; rl->view_dir = pole; int i; for(i = 0; i < num_comp; i++){ r3_t v = normals[i]; rl->bradius[i] = span/(resolution*v.c[2]); rl->coeficients[i] = 1.0; // fprintf(stderr, " %3d ", i); //r3_print(stderr, &((*bcenter)[i])); //fprintf(stderr, " %9.6f\n", br[i]); } free(pontos); return rl; } double radialbasis_phi(int r, double* x, void* l_data){ r3_t normal = (r3_t) {{x[0],x[1],x[2]}}; radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; r3_t bcenter = rl->bcenter[r]; double bradius = rl->bradius[r]; double dot = r3_dot(&normal,&bcenter); double d2 = (1.0 - dot)/(bradius*bradius); if(d2 >= 1.0) return 0.0; return (1 - d2)*(1 - d2); } int radialbasis_get_number_components(void* l_data){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; return rl->num_comp; } void radialbasis_set_alphas(double* C, void* l_data){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; for(i = 0; i < rl->num_comp; i++){ rl->coeficients[i] = C[i]; } } void radialbasis_get_alphas(void* l_data,double* C,int n){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; assert(rl->num_comp == n); int i; for(i = 0; i < rl->num_comp; i++){ C[i] = rl->coeficients[i]; } } void* radialbasis_copy_lighting_data(void* l_data){ radialbasis_lighting_t* pl = (radialbasis_lighting_t*) l_data; radialbasis_lighting_t* npl = (radialbasis_lighting_t*) malloc(sizeof(radialbasis_lighting_t)); npl->num_comp = pl->num_comp; npl->bcenter= (r3_t*)malloc(sizeof(r3_t)*(npl->num_comp)); npl->coeficients = (double*)malloc(sizeof(double)*(npl->num_comp)); npl->bradius = (double*)malloc(sizeof(double)*(npl->num_comp)); int i; for(i = 0; i < npl->num_comp;i++){ npl->bradius[i] = pl->bradius[i]; npl->bcenter[i] = pl->bcenter[i]; npl->coeficients[i] = pl->coeficients[i]; } return npl; } void radialbasis_release_lighting_data(void* l_data){ radialbasis_lighting_t* pl = (radialbasis_lighting_t*) l_data; free(pl->bradius); free(pl->bcenter); free(pl->coeficients); free(l_data); } double radialbasis_shading(double* x,void* l_data){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; int r; double val = 0; for(r = 0; r < rl->num_comp; r++){ val+= (rl->coeficients[r]*radialbasis_phi(r,x,l_data)); } return val; } void radialbasis_write_params(FILE* arq,void* l_data){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data; fprintf(arq," %d\n",rl->num_comp); int i; for(i = 0; i < rl->num_comp;i++){ fprintf(arq,"%8.6lf %8.6lf",rl->coeficients[i],rl->bradius[i]); r3_gen_print(arq, &(rl->bcenter[i]), "%9.6lf", " ", " ", ""); fprintf(arq,"\n"); } } void* radialbasis_read_params(FILE* arq){ radialbasis_lighting_t* rl = (radialbasis_lighting_t*)malloc(sizeof(radialbasis_lighting_t)); fscanf(arq,"%d",&(rl->num_comp)); int i; int num_comp = rl->num_comp; rl->bradius = (double*)malloc(sizeof(double)*num_comp); rl->coeficients = (double*)malloc(sizeof(double)*num_comp); rl->bcenter = (r3_t*)malloc(sizeof(r3_t)*num_comp); for(i = 0; i < rl->num_comp;i++){ fscanf(arq,"%lf %lf",&(rl->coeficients[i]),&(rl->bradius[i])); //r3_gen_print(arq, &(rl->bcenter[i]), "%9.6lf", " ", " ", ""); fscanf(arq,"%lf %lf %lf",&(rl->bcenter[i].c[0]),&(rl->bcenter[i].c[1]),&(rl->bcenter[i].c[2])); // fprintf(arq,"\n"); } return rl; } void radialbasis_update_weights(void* l_data,double** X, double* F,double* weights, int n,double * sigma){ double sumd2 = 0.0; double sumw = 1.0e-200; double sigma2 = *sigma; sigma2 = sigma2*sigma2; int i; for(i = 0; i< n; i++){ // r3_t normal = (r3_t){{X[i][0],X[i][1],X[i][2]}}; double computedGO = radialbasis_shading(X[i],l_data); double originalGO = F[i]; double d2 = (computedGO-originalGO)*(computedGO-originalGO); weights[i] = (d2 > 5.0*sigma2 ? 0.0 : exp(-d2/(2.0*sigma2)) ); sumd2+= d2; sumw+=weights[i]; } *sigma = sqrt(sumd2/(double)n); } double radialbasis_compare_lightings(void* l_data1,void* l_data2){ return 0.0; } /*Stereographic projection of normals based on polynoms*/ harmonicSG_lighting_t* harmonicSG_init_components(int g,r3_t view_dir){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*)malloc(sizeof(harmonicSG_lighting_t)); int num_comp = (g+1)*(g+2)/2.0; pl->num_comp = num_comp; pl->view_dir = view_dir; pl->coeficients = (double*)malloc(sizeof(double)*num_comp); pl->tuples = (double**)malloc(sizeof(double*)*num_comp); int k,i,j; for(i = 0; i < num_comp; i++){ pl->coeficients[i] = 1.0; pl->tuples[i] = (double*)malloc(sizeof(double)*2); } int count = 0; for(k = 0; k <= g; k++){ for(i = k; i >= 0; i--){ j = k - i ; pl->tuples[count][0] = i; pl->tuples[count][1] = j; count++; } } pl->num_comp = count; return pl; } double harmonicSG_phi(int r, double* x, void* l_data){ /*Stereographic coordinates*/ double X,Y; X = x[0]/(x[2] + 1.0); Y = x[1]/(x[2] + 1.0); harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; int ri = pl->tuples[r][0]; int rj = pl->tuples[r][1]; double phiR = (pow(X,ri))*(pow(Y,rj)); //phiR = phiR*pl->weights[r]; return phiR; } int harmonicSG_get_number_components(void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; return pl->num_comp; } void harmonicSG_set_alphas(double* C,void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; for(i = 0; i < pl->num_comp; i++){ pl->coeficients[i] = C[i]; } } double harmonicSG_get_alphas(void* l_data,double* C,int n){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; assert(pl->num_comp == n); int i; for(i = 0; i < pl->num_comp; i++){ C[i] = pl->coeficients[i]; } } void* harmonicSG_copy_lighting_data(void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; harmonicSG_lighting_t* npl = (harmonicSG_lighting_t*) malloc(sizeof(harmonicSG_lighting_t)); npl->num_comp = pl->num_comp; npl->tuples = (double**)malloc(sizeof(double*)*(npl->num_comp)); npl->coeficients = (double*)malloc(sizeof(double)*(npl->num_comp)); int i; for(i = 0; i < npl->num_comp;i++){ npl->coeficients[i] = pl->coeficients[i]; npl->tuples[i] = (double*)malloc(sizeof(double)*2); int j; for(j = 0; j < 2;j++){ npl->tuples[i][j] = pl->tuples[i][j]; } } return npl; } void harmonicSG_release_lighting_data(void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; int i; for(i = 0; i < pl->num_comp;i++){ free(pl->tuples[i]); } free(pl->tuples); free(pl->coeficients); free(l_data); } double harmonicSG_shading(double* x,void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; int r; double val = 0; for(r = 0; r < pl->num_comp; r++){ val+= (pl->coeficients[r]*harmonicSG_phi(r,x,l_data)); } return val; } void harmonicSG_write_params(FILE* arq,void* l_data){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) l_data; fprintf(arq," %d\n",pl->num_comp); int i; for(i = 0; i < pl->num_comp;i++){ fprintf(arq, "%8.6lf %8.6lf %8.6lf\n",pl->coeficients[i],pl->tuples[i][0],pl->tuples[i][1] ); } } void* harmonicSG_read_params(FILE* arq){ harmonicSG_lighting_t* pl = (harmonicSG_lighting_t*) malloc(sizeof(harmonicSG_lighting_t)); fscanf(arq,"%d",&(pl->num_comp)); int i; int num_comp = pl->num_comp; pl->coeficients = (double*)malloc(sizeof(double)*num_comp); pl->tuples = (double**)malloc(sizeof(double*)*num_comp); for(i = 0; i < pl->num_comp;i++){ pl->tuples[i] = (double*)malloc(sizeof(double)*2); fscanf(arq, "%lf %lf %lf",&(pl->coeficients[i]),&(pl->tuples[i][0]),&(pl->tuples[i][1])); } return pl; } void harmonicSG_update_weights(void* l_data,double** X, double* F,double* weights, int n,double * sigma){ double sumd2 = 0.0; double sumw = 1.0e-200; double sigma2 = *sigma; sigma2 = sigma2*sigma2; int i; for(i = 0; i< n; i++){ // r3_t normal = (r3_t){{X[i][0],X[i][1],X[i][2]}}; double computedGO = harmonicSG_shading(X[i],l_data); double originalGO = F[i]; double d2 = (computedGO-originalGO)*(computedGO-originalGO); weights[i] = (d2 > 5.0*sigma2 ? 0.0 : exp(-d2/(2.0*sigma2)) ); sumd2+= d2; sumw+=weights[i]; } *sigma = sqrt(sumd2/(double)n); } double harmonicSG_compare_lightings(void* l_data1,void* l_data2){ return 0.0; } /*Compact Light model*/ compact_lighting_t* compact_init_components( r3_t dirlight, double err_dirligh, double radlight, double err_radlight, double shine, double err_shine, bool_t has_isotropic, r3_t backnormal, double err_backnormal, r3_t view_dir ){ compact_lighting_t* cl = (compact_lighting_t*)malloc(sizeof(compact_lighting_t)); cl->dirlight = dirlight; cl->err_dirligh = err_dirligh; if( r3_norm(&dirlight) == 0){ cl->Hc = FALSE; cl->dirlight.c[0] = 0; cl->dirlight.c[1] = 0; cl->dirlight.c[2] = 1; }else{ cl->Hc = TRUE; } cl->radlight = radlight; cl->err_radlight = err_radlight; cl->shine = shine; cl->err_shine = err_shine; if( shine == INFINITY){ cl->Hl = FALSE; }else{ cl->Hl = TRUE; } cl->backnormal = backnormal; if( r3_norm(&backnormal) == 0){ cl->Hf = FALSE; cl->backnormal.c[0] = 0; cl->backnormal.c[1] = 0; cl->backnormal.c[2] = 1; }else{ cl->Hf = TRUE; } cl->Hi = has_isotropic; cl->view_dir = view_dir; cl->Ec = cl->El = cl->Ei = cl->Ef = 0; return cl; } void compact_determine_basis_indices(void* l_data,int* ne, int* num_compac, int* num_phong, int* num_isotropic, int* num_backplane){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int n = 0; if((num_compac!= NULL) && cl->Hc){ *num_compac = n; n++; }else{num_compac = -1;} if((num_phong!= NULL) && cl->Hl){ *num_phong = n; n++ }else{ num_phong = -1 ;} if((num_isotropic!= NULL) && cl->Hi){ *num_isotropic= n; n++ }else{ num_isotropic = -1 ;} if((num_backplane!= NULL) && cl->Hl){ *num_backplane = n; n++ }else{ num_backplane = -1 ;} *ne = n; } r3_t compute_bissetriz(r3_t u, r3_t v); r3_t compute_bissetriz(r3_t u, r3_t v){ r3_t w; /*w = 2(v.u)v - u*/ r3_scale(2*r3_dot(&v,&u),&v,&w); r3_sub(&w,&u,&w); return w; } double compact_phi(int r, double* x,void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; r3_t normal = (r3_t){{x[0],x[1],x[2]}}; auto double Phong(void); double Phong(void){ r3_t w = compute_bissetriz(&(cl->dirlight),&normal); double phong_base = fmax(0, r3_dot(&w,&(cl->view_dir))); return pow(phong_base,cl->shine) ; } int n_comp; int num_compac,num_phong,num_isotropic,num_backplane; compact_determine_basis_indices(l_data,&n_comp,&num_compac,&num_phong,&num_isotropic,&num_backplane); if(num_compac == r){ return pst_lamp_geom_factor(&normal, &(cl->dirlight),cl->radlight); }else if (num_phong == r){ return Phong(); }else if (num_isotropic == r){ return 1; }else if( num_backplane == r){ return 1 - r3_dot(&normal,&(cl->backnormal)); } fprintf(stderr,"libraab - compact_phi: ERROR - invalid value for r = %d within range [0-%d]\n",r,n_comp -1 ); assert(FALSE); } int compact_get_number_components(void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int n_comp; compact_determine_basis_indices(l_data,&n_comp,NULL,NULL,NULL,NULL); return ne; } void compact_set_alphas(double* C,int n, void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int n_comp; int num_compac,num_phong,num_isotropic,num_backplane; compact_determine_basis_indices(l_data,&n_comp,&num_compac,&num_phong,&num_isotropic,&num_backplane); assert(n_comp == n); cl->Ec = (num_compac >= 0 ? C[num_compac] : 0); cl->El = (num_phong >= 0 ? C[num_phong] : 0); cl->Ei = (num_isotropic >= 0 ? C[num_isotropic] : 0); cl->Ef = (num_backplane >= 0 ? C[num_backplane] : 0); } void compact_get_alphas(void* l_data,double* C,int n){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int n_comp; int num_compac,num_phong,num_isotropic,num_backplane; compact_determine_basis_indices(l_data,&n_comp,&num_compac,&num_phong,&num_isotropic,&num_backplane); assert(n_comp == n); if( num_compac >= 0 ) C[num_compac] = cl->Ec; if( num_phong >= 0 ) C[num_phong] = cl->El; if( num_isotropic >= 0 ) C[num_isotropic] = cl->Ei; if( num_backplane >= 0 ) C[num_backplane] = cl->Ef; } void* compact_copy_lighting_data(void* l_data){ compact_lighting_t* cl1 = (compact_lighting_t*)l_data; compact_lighting_t* cl2 = (compact_lighting_t*) malloc(sizeof(compact_lighting_t)); *cl2 = *cl1; return cl2; } void compact_release_lighting_data(void* l_data){ free(l_data); } double compact_shading(double* x,void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int n_comp = compact_get_number_components(l_data); int i; double val = 0; for(i = 0; i > n_comp; i++){ val+= compact_phi(x,i,l_data); } return val ; } void compact_write_params(FILE* arq,void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; fprintf(arq,"%9.6lf %9.6lf %9.6lf %9.6lf %9.6lf %9.6lf %9.6lf ",cl->Ec,cl->rho,cl->dirlight.c[0],cl->dirlight.c[1],cl->dirlight.c[2],cl->radlight, cl->err_dirligh ); fprintf(arq,"%9.6lf %9.6lf %9.6lf ",cl->El,cl->shine,cl->err_shine); fprintf(arq,"%9.6lf ",cl->Ei); fprintf(arq,"%9.6lf %9.6lf %9.6lf %9.6lf %9.6lf ",cl->Ef,cl->backnormal.c[0],cl->backnormal.c[1],cl->backnormal.c[2],cl->err_backnormal); fprintf(arq,"%d %d %d %d", cl->Hc, cl->Hl, cl->Hi, cl->Hf); fprintf(arq,"\n"); } void* compact_read_params(FILE* arq){ compact_lighting_t* cl = (compact_lighting_t*)malloc(sizeof(compact_lighting_t)); fscanf(arq,"%lf %lf %lf %lf %lf %lf %lf",&(cl->Ec),&(cl->rho),&(cl->dirlight.c[0]),&(cl->dirlight.c[1]),&(cl->dirlight.c[2]),&(cl->radlight),&(cl->err_dirligh) ); fscanf(arq,"%lf %lf %lf",&(cl->El),&(cl->shine),&(cl->err_shine)); fscanf(arq,"%lf",&(cl->Ei)); fscanf(arq,"%lf %lf %lf %lf %lf",&(cl->Ef),&(cl->backnormal.c[0]),&(cl->backnormal.c[1]),&(cl->backnormal.c[2]),&(cl->err_backnormal)); fscanf(arq,"%d %d %d %d",&(cl->Hc),&(cl->Hl),&(cl->Hi),&(cl->Hf)); return cl; } void compact_update_weights(void* l_data,double** X, double* F,double* weights, int n,double* sigma){ compact_lighting_t* cl = (compact_lighting_t*)l_data; int i; r3_t u = cl->dirlight; double K = cl->shine; double rho = cl->radlight; double eps = cl->err_dirligh; for(i = 0; i < n; i++){ r3_t normal = (r3_t){{X[i][0],X[i][1],X[i][2]}}; /*Verifies if it is a region of self-shadow or terminator*/ double wt = (r3_dot(&u,&normal) < cos(rho + eps) ? 0 : 1); /*Check if it is in region of phong highlight*/ double w_bissetriz = compute_bissetriz(u,normal); double tol = 10e-5; double lambda = acos(pow(tol,1/K)); double wl = ( r3_dot(&normal,&w_bissetriz) > cos( (eps/2.0) + lambda) ? 0 : 1); weights[i] = wl*wt; } } double compact_compare_lightings(void* l_data1,void* l_data2){ compact_lighting_t* cl1 = (compact_lighting_t*)l_data1; compact_lighting_t* cl2 = (compact_lighting_t*)l_data2; if(cl1->Hc != cl2->Hc){ return INFINITY;} if(cl1->Hl != cl2->Hl){ return INFINITY;} if(cl1->Hi != cl2->Hi){ return INFINITY;} if(cl1->Hf != cl2->Hf){ return INFINITY;} int n_comp1 = compact_get_number_components(l_data1); double v1[n_comp1]; compact_get_alphas(l_data1,v1,n_comp1); int n_comp2 = compact_get_number_components(l_data2); double v2[n_comp2]; compact_get_alphas(l_data2,v2,n_comp2); assert(n_comp1 == n_comp2); return rn_dist(n_comp1,v1,v2); } int compact_num_nl_params(void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)malloc(sizeof(compact_lighting_t)); /*Only valid params*/ num_parms = 0; if(cl->Hc) num_parms+= 3; if(cl->Hl) num_parms+= 1; if(cl->Hf) num_parms+= 2; return num_parms; } void convert_r3_to_polar(r3_t v, double* rho_v, double* theta_v,r3_t view_dir); void convert_r3_to_polar(r3_t v, double* rho_v, double* theta_v,r3_t view_dir){ r3_t v_novo; r3x3_t VM = compute_normal_correction_matrix(view_dir); r3x3_map_row(&v, &VM, &v_novo); double X = v_novo.c[0]/(1 + v_novo.c[2]); double Y = v_novo.c[1]/(1 + v_novo.c[2]); *rho_v = sqrt( X*X + Y*Y); *theta_v = atan2(Y,X); } r3_t convert_polar_to_r3(double rho, double theta, r3_t view_dir); r3_t convert_polar_to_r3(double rho, double theta, r3_t view_dir){ r3_t v; double X = rho*cos(theta); double Y = rho*sin(theta); r3_t v_novo; double den = (1+ X*X + Y*Y); v_novo.c[0] = 2*X/den; v_novo.c[1] = 2*Y/den; v_novo.c[2] = (1 - X*X - Y*Y)/den; r3x3_t VM = compute_normal_correction_matrix(view_dir); r3x3_inv(&VM,&VM); r3x3_map_row(&v_novo, &VM, &v); return v; } void compact_pack_parameter(void* l_data,double* parameters,int n){ int num_parms = compact_num_nl_params(l_data); compact_lighting_t* cl = (compact_lighting_t*)l_data; r3_t view_dir = cl->view_dir; int ind_parm = 0; /*primeiro o vetor de direção de luz*/ if(cl->Hc){ r3_t u = cl->dirlight; double rho_u,theta_u; convert_r3_to_polar(u,&rho_u,&theta_u); parameters[ind_parm] = rho_u; parameters[ind_parm+1] = theta_u; parameters[ind_parm+2] = cl->radlight; ind_parm+=3; } /*Now K and rho*/ if(cl->Hl){ parameters[ind_parm] = cl->shine; ind_parm+=1; } /*Now Xi*/ if( cl->Hf){ r3_t xi = cl->backnormal; double rho_xi,theta_xi; convert_r3_to_polar(xi,&rho_xi,&theta_xi); parameters[ind_parm] = rho_xi; parameters[ind_parm+1] = theta_xi; ind_parm+=2; } int i; for(i = 0; i < num_parms ; i++){ assert(!isnan(parameters[i])); } } void compact_unpack_parameter(double* parameters,int n,void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; num_parms = compact_num_nl_params(l_data); assert(n == num_parms); int i; for(i = 0; i < num_parms; i++){ assert(!isnan(parameters[i])); } /*Copy unpacked data*/ int ind_parm = 0; if(cl->Hc){ double rho_u = parameters[ind_parm]; double theta_u = parameters[ind_parm+1]; cl->dirlight = convert_polar_to_r3(rho_u,theta_u); cl->radlight = parameters[ind_parm+2]; ind_parm+=3; } if( cl->Hl){ cl->shine = parameters[ind_parm]; ind_parm+=1; } if( cl->Hf){ double rho_xi = parameters[ind_parm]; double theta_xi = parameters[ind_parm+1]; cl->backnormal = convert_polar_to_r3(rho_xi,theta_xi); ind_parm+=2; } return; } r3_t find_ortho_vector(r3_t view_dir,r3_t u); /*Find a vector orthogonal {t} to view_dir and coplanar with u*/ r3_t find_ortho_vector(r3_t view_dir,r3_t u){ r3_t t,w; r3_cross(&view_dir,&u,&w); double nw = r3_dir(&w,&w); if(nw == 0){ t = (r3_t){{1,0,0}}; }else{ r3_cross(&view_dir,&w,&t); double nt = r3_dir(&t,&t); } return t; } r3_t bend_towards(r3_t v, r3_t u, double ang); /*Bends a vector {u} towards unit vector {v} by {ang} radians */ r3_t bend_towards(r3_t v, r3_t u, double ang){ r3_t t = find_ortho_vector(view_dir,u); double uv = r3_dot(&v,&u); double ut = r3_dot(&t,&u); double u1v = +uv*cos(ang) + ut*sin(ang); double u1t = -uv*sin(ang) + ut*cos(ang); r3_t r; r3_mix(u1v,&v,u1t,&t,&r); return r; } void convert_direction_and_error_to_polar(r3_t u, r3_t view_dir, double err_u, double* rho_u, double* theta_u, double* err_rho_u, double* err_theta_u); void convert_direction_and_error_to_polar(r3_t u, r3_t view_dir, double err_u, double* rho_u, double* theta_u, double* err_rho_u, double* err_theta_u){ } double compact_pack_nl_param_errors(void* l_data,double* par_errors, int n){ compact_lighting_t* cl = (compact_lighting_t*)l_data; r3_t view_dir = cl->view_dir; int ind_parm = 0; /*primeiro o vetor de direção de luz*/ if(cl->Hc){ double err_rho_u = 0; double err_theta_u = 0; if(cl->err_dirligh > 0) { r3_t u = cl->dirlight; double rho_u,theta_u; convert_r3_to_polar(u,&rho_u, &theta_u,view_dir); /*Find the two vectors {u1,u2} towards and away from ${view_dir} by angle {err_dirligh}*/ r3_t u1 = bend_towards(view_dir,u,+cl->err_dirligh); r3_t u2 = bend_towards(view_dir,u,-cl->err_dirligh); double rho_u1,theta_u1; convert_r3_to_polar(u1,&rho_u1, &theta_u1,view_dir); double rho_u2,theta_u2; convert_r3_to_polar(u2,&rho_u2, &theta_u2,view_dir); if( fabs(theta_u1 - theta_u2) < 1.0e-5){ err_rho_u = fabs(rho_u1 - rho_u2)/2.0; err_theta_u = atan2(err_rho_u,sqrt(fmax(0,rho_u*rho_u - err_rho_u*err_rho_u)) ); }else if (fabs( fabs(theta_u1 - theta_u2) - M_PI) < 10e-5 ) { err_rho_u = fabs(rho_u1 + rho_u2)/2.0; err_theta_u = M_PI; }else{ /*Assume that one {u1,u2} is very close to {view_dir}*/ erro_rho_u = fmax(rho_u1,rho_u2); err_theta_u = M_PI/2.0; } } parameters[ind_parm] = err_rho_u; parameters[ind_parm+1] = err_theta_u; parameters[ind_parm+2] = cl->err_radlight; ind_parm+=3; } /*Now K and rho*/ if(cl->Hl){ parameters[ind_parm] = cl->shine; ind_parm+=1; } /*Now Xi*/ if( cl->Hf){ r3_t xi = cl->backnormal; double rho_xi,theta_xi; convert_r3_to_polar(xi,&rho_xi,&theta_xi); parameters[ind_parm] = rho_xi; parameters[ind_parm+1] = theta_xi; ind_parm+=2; } int i; for(i = 0; i < num_parms ; i++){ assert(!isnan(parameters[i])); } } double compact_unpack_nl_param_errors(double* par_errors,int n, void* l_data){ compact_lighting_t* cl = (compact_lighting_t*)l_data; compact_lighting_t* cl = (compact_lighting_t*)l_data; nl_parms = compact_num_nl_params(l_data); assert(nl_parms == n); }