#include <models.h>
#include <assert.h>
#include <tabela.h>

/*Generic model type*/

lighting_model_t* create_lighting_model(lighting_type_t type){
  
  lighting_model_t* lm = (lighting_model_t*) malloc(sizeof(lighting_model_t));
  lm->lighting_data = NULL;
  lm->type = type;
  if(type == POINTLIKE_MODEL){
    lm->phi = pointlike_phi;
    lm->get_num_components = pointlike_get_number_components;
    lm->retrieve_components = pointlike_retrieve_components;
    lm->update_weights = pointlike_update_weights;
    lm->validate_results= NoValidation;
    lm->shading = pointlike_shading;
    lm->write_param = pointlike_write_params;
    lm->copy_lighting_data =  pointlike_copy_lighting_data;
    lm->compare_lightings = pointlike_compare_lightings;
    lm->release_lighting_data = pointlike_release_lighting_data;
  }
  else if(type == HARMONIC_MODEL){
    lm->phi = harmonic_phi;
    lm->get_num_components = harmonic_get_number_components;
    lm->retrieve_components = harmonic_retrieve_components;
    lm->update_weights = harmonic_update_weights;
    lm->validate_results= NoValidation;
    lm->shading = harmonic_shading;
    lm->write_param = harmonic_write_params;
    lm->copy_lighting_data =  harmonic_copy_lighting_data;
    lm->compare_lightings = harmonic_compare_lightings;
    lm->release_lighting_data = harmonic_release_lighting_data;
  }else if(type == RADIALBASIS_MODEL){
    lm->phi = radialbasis_phi;
    lm->get_num_components = radialbasis_get_number_components;
    lm->retrieve_components = radialbasis_retrieve_components;
    lm->update_weights = radialbasis_update_weights;
    lm->validate_results= NoValidation;
    lm->shading = radialbasis_shading;
    lm->write_param = radialbasis_write_params;
    lm->copy_lighting_data =  radialbasis_copy_lighting_data;
    lm->compare_lightings = radialbasis_compare_lightings;
    lm->release_lighting_data = radialbasis_release_lighting_data;
  }else{
    free(lm);
    return NULL;
  }
  return lm;
}

/*Lambert Model*/

double pointlike_phi(int r, r3_t *normal,void* l_data){
  
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  if(r < 3){ /* Elements for direct lighting by point source: */
    double s = r3_dot(&(ll->f), normal);
    return (s <= 0 ? 0 : normal->c[r]);
  }else if (r == 3){ /* Element for indirect lighing: */
    return (ll->BPL ? 0.5*(1 - normal->c[2]) : 1);
  }else{ /* No such element: */
    assert(FALSE);
  }
  
}

int pointlike_get_number_components(void* l_data){
  return 4;
}

pointlike_lighting_t* pointlike_init_components(double albedo,bool_t BPL){
  
  pointlike_lighting_t* ll = (pointlike_lighting_t*) malloc(sizeof(pointlike_lighting_t));
  ll->f =  (r3_t){{ 0,0,1}};
  ll->Ed = 0;
  ll->Eb = 0;
  ll->Ei = 0;
  ll->BPL = BPL;
  ll->albedo = albedo;
  return ll;
  
}

void pointlike_retrieve_components(double* C, int n,void* l_data){
  
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  ll->f.c[0] = C[0];
  ll->f.c[1] = C[1];
  ll->f.c[2] = C[2];
  ll->Ed = r3_dir(&(ll->f),&(ll->f))/ll->albedo;
  if(ll->BPL){
    ll->Eb = C[3]/ll->albedo;
    ll->Ei = 0;
  }else{
    ll->Eb = 0;
    ll->Ei = C[3]/ll->albedo;
  }
  
}

void pointlike_update_weights(void* l_data,r3_t* X, double* F,double* weights, int n,double* sigma){
  int i;
  
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  for(i = 0; i < n; i++){
    r3_t normal = X[i];
    double s = r3_dot(&(ll->f), &normal);
    if(s > 0){
      weights[i] = 1.0;
    }else{
      weights[i] = 0.0;
    }
  }
  
}


double pointlike_shading(r3_t* normal,void* l_data){
  
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  double vald = fmax(0,ll->Ed*r3_dot(&(ll->f),normal));
  double valb = ll->Eb*(1-normal->c[3])/2.0;
  double vali = ll->Ei;
  double radiation = ll->albedo*(vald+valb+vali);
  return radiation;
  
}

void pointlike_write_params(FILE* arq,void* l_data){
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  fprintf(arq, "  %8.6lf", ll->albedo);
  r3_gen_print(arq, &(ll->f), "%9.6lf", "  ", " ", "");
  fprintf(arq, "  %8.6lf %8.6lf %8.6lf", ll->Ed, ll->Eb, ll->Ei);
  fprintf(arq, "\n");
}


void* pointlike_copy_lighting_data(void* l_data){
  pointlike_lighting_t* ll = (pointlike_lighting_t*) l_data;
  pointlike_lighting_t* npl = (pointlike_lighting_t*) malloc(sizeof(pointlike_lighting_t));
  npl->f = ll->f;
  npl->albedo = ll->albedo;
  npl->Ei = ll->Ei;
  npl->Ed = ll->Ed;
  npl->Eb = ll->Eb;
  npl->BPL = ll->BPL;
  return npl;
}

double pointlike_compare_lightings(void* l_data1,void* l_data2){
  pointlike_lighting_t* ll1 = (pointlike_lighting_t*) l_data1;
  pointlike_lighting_t* ll2 = (pointlike_lighting_t*) l_data2;
  r3_t f = ll1->f;
  r3_t f_new = ll2->f;
  double diff = r3_dist(&f, &f_new);
  return diff;
}

void pointlike_release_lighting_data(void* l_data){
  free(l_data);
}

/*Polynomial model*/


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] );
  }
  
}

double harmonic_phi(int r, r3_t* normal, void* l_data){
  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;
}


harmonic_lighting_t* harmonic_init_components(int g){
  
  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->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;
  
}

void harmonic_retrieve_components(double* C, int n,void* l_data){
  harmonic_lighting_t* pl = (harmonic_lighting_t*) l_data;
  int i;
  assert(pl->num_comp >= n);
  for(i = 0; i < n; i++){
    pl->coeficients[i] = C[i];
  }
  for(i = i ; i < pl->num_comp; i++){
    pl->coeficients[i] = 0;
  }
  
}



double harmonic_shading(r3_t* normal,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,normal,l_data));
  }
  
  return val;
}


void harmonic_update_weights(void* l_data,r3_t* 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 = X[i];
    double computedGO = harmonic_shading(&normal,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);
}

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));
  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;
}

double harmonic_compare_lightings(void* l_data1,void* l_data2){
  return 0.0;
}

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);
}


/*Radial Basis Model*/



double radialbasis_phi(int r, r3_t* normal, void* l_data){
  
  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);
  
}

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");
  }
  
}

int radialbasis_get_number_components(void* l_data){
  radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data;
  return rl->num_comp;
}


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;

	
	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;
}

void radialbasis_retrieve_components(double* C, int n,void* l_data){
  
  radialbasis_lighting_t* rl = (radialbasis_lighting_t*) l_data;
  int i;
  assert(rl->num_comp >= n);
  for(i = 0; i < n; i++){
    rl->coeficients[i] = C[i];
  }
  for(i = i ; i < rl->num_comp; i++){
    rl->coeficients[i] = 0;
  }
  
}

void radialbasis_update_weights(void* l_data,r3_t* 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 = X[i];
    double computedGO = radialbasis_shading(&normal,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_shading(r3_t* normal,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,normal,l_data));
  }
  
  return val;
}

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;
}

double radialbasis_compare_lightings(void* l_data1,void* l_data2){
  return 0.0;
}

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);
}


/*Validation functions*/
int NoValidation(double* c, int basis_size, bool_t* validity_array){
  return basis_size;
}