/* Programa que faz o casamento do modelo com os dados obtidos.
   O arquivo de entrada dos vertices deve estar no formato
   x y z    para cada vertice.  */

#include <stdio.h>
#include <math.h>

#define NUM_CASA 3 /* num minimo de vertices casados (dados com modelo transladado) */

#define SIZEX 512
#define SIZEY 512

float baseline, foco, const_cord_x, const_cord_y;
float deltamax;
float **mv,  /* matriz de vertices do modelo */
      **mu_3d,  /* matriz de vertices dos dados em 3d */
      **mu_iso; /* matriz de vertices dos dados isotropicos */
float **mv_homog;  /* mat. para coord. homogeneas dos
			   vertices do modelo*/
float **mv_map_hom,  /* vertices do modelo mapeados pela 
	                          mat final em coord. homogeneas */
      **mv_map_iso; /* vertices mv_map_hom em 
	                    cord. cartezianas isotropicos */


FILE  *saida, *file_log;

int nvertice_dados, nvertice_mod;  /* num. de vertices obtidos dos dados
                                      e do modelo */


/* multiplicacao de vetor linha por matriz  em R4 */
void prod_row_mat(x, m, res)
float *x, m[4][4], *res; /* x[4], res[4] */
{
    float s;
    int j, k;

    for(j=0; j<4 ;j++)
      {
	  s = 0.0;
	  for(k=0; k<4 ;k++)
	    {
		s += x[k]*m[k][j];
	    }
	  res[j] = s;
      }
}


/* determinante de uma matriz 3x3 */
float det3x3(m)
float m[3][3];
{
    float res, d0, d1, d2;

    d0 = m[1][1]*m[2][2] - m[1][2]*m[2][1];
    d1 = m[1][0]*m[2][2] - m[1][2]*m[2][0];
    d2 = m[1][0]*m[2][1] - m[1][1]*m[2][0];
    res = d0 * m[0][0] - d1 * m[0][1] + d2 * m[0][2];

    return(res);
}

/* cofator do elemento [ix][jx] de uma matriz 4x4 */
float cof(matin, ix, jx)
float matin[4][4];
int ix, jx;
{
  float t[3][3];
  float res, d;
  int ii, jj, i, j;

  ii = 0;
  for(i=0; i<4; i++)
    {
	if(i!=ix)
	  {
	      jj = 0;
	      for(j=0; j<4; j++)
		{
		    if(j!=jx)
		      {
			  t[ii][jj] = matin[i][j];
			  jj++;
		      }
		}
	      ii++;
	  }
    }

  d = det3x3(t);

  if(((ix+jx)%2)!=0)
    res = -d;
  else
    res = d;

  return(res);
}


/* determinante de uma matriz 4x4 */
float det4x4(m)
float m[4][4];
{
  float res;

  res = m[0][0] * cof( m, 0, 0 ) + m[0][1] * cof( m, 0, 1 ) + 
        m[0][2] * cof( m, 0, 2 ) + m[0][3] * cof( m, 0, 3 );

  return(res);
}
		  


/* inversa de uma matriz 4x4 : m * result = I * det(m)**2   */
void inversa(matin, matout)
float matin[4][4], matout[4][4];
{
  float d;
  int i, j;

  for(i=0; i<4; i++)
    for(j=0; j<4; j++)
      matout[i][j] = 0.0;
  
  d = det4x4(matin);
  if(d == 0.0)
    return;

  for(i=0; i<4; i++)
    for(j=0; j<4; j++)
      matout[i][j] = cof(matin, j, i)/d;

  return;
}




void mat_prod(a, b, o)  /* produto de matrizes o = ab
				       a(4 x 4),  b(4 x 4) */
float a[4][4], b[4][4], o[4][4];
{
    int i, j, k;
    float s;

    for(k=0; k<4; k++)
      for(i=0; i<4; i++)
	{
	    s = 0.0;
	    for(j=0; j<4; j++)
	      {
		  s += a[k][j]*b[j][i];
	      }
	    o[k][i] = s;
	}
}
	
	
 
     

void crossp(a, b, o) /* produto vetorial o = a x b  de vetores R3  */
float *a, *b, *o;
{
    float d;
    o[0] = a[1]*b[2] - a[2]*b[1];
    o[1] = a[2]*b[0] - a[0]*b[2];
    o[2] = a[0]*b[1] - a[1]*b[0];

    d = sqrt(o[0]*o[0] + o[1]*o[1] + o[2]*o[2]);
    o[0] /= d;
    o[1] /= d;
    o[2] /= d;
}

float normalize(a) /* calcula o modulo do vetor em R3  */
float *a;
{
    float d;

    d = sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);

    a[0] /= d;
    a[1] /= d;
    a[2] /= d;

    return(d);
}

/* dist() - calcula a distancia entre dois ptos no espaco R3.
            Retorna a distancia calculada.   */
float dist(p1, p2)
float *p1,  /* ponto 1 ; vetor[3] */
      *p2;  /* ponto 2 ; vetor[3] */
{
 int i;
 float dist = 0.0;  /* distancia entre p1 e p2  */

 for(i=0; i<3; i++)
   {
     dist += (p1[i] - p2[i]) * (p1[i] - p2[i]);
   }

 dist = sqrt(dist);

 return(dist);
}

/* para visualizar uma matriz */
void ver2(m, nlinhas, ncol)
int nlinhas, ncol;
float m[nlinhas][ncol];
{
  int i, j;
  printf(stderr,"\n");
  for(i = 0; i<nlinhas; i++)
    {
      for(j=0; j<ncol; j++)
	{
	  fprintf(file_log,"   %f",m[i][j]);
        }

      fprintf(file_log,"\n\n");
    }
}


/* para visualizar uma matriz- forma diferente de alocacao*/
void ver1(m, nlinhas, ncol)
int nlinhas, ncol;
float **m;
{
  int i, j;
  printf(stderr,"\n");
  for(i = 0; i<nlinhas; i++)
    {
      for(j=0; j<ncol; j++)
	{
	  fprintf(file_log,"   %f",m[i][j]);
        }

      fprintf(file_log,"\n\n");
    }
}

/* funcao conversao de iso para 3d */
void conv_iso_3d(x_iso, y_iso, z_iso, x_3d, y_3d, z_3d)
float x_iso, y_iso, z_iso;
float *x_3d, *y_3d, *z_3d;
{
  *x_3d = baseline * (x_iso - const_cord_x)/(2.0*z_iso);
  *y_3d = baseline * (y_iso*1.414214 - const_cord_y)/(2.0*z_iso);
  *z_3d = foco*baseline/(2.0*z_iso);
}  

/* funcao conversao de 3d para iso */
void conv_3d_iso(x_3d, y_3d, z_3d, x_iso, y_iso, z_iso)
float *x_iso, *y_iso, *z_iso;
float x_3d, y_3d, z_3d;
{
  *z_iso = foco*baseline/(2.0*z_3d);
  *x_iso = 2.0*(*z_iso)*x_3d/baseline + const_cord_x;
  *y_iso = 1.414214*(*z_iso)*y_3d/baseline + 0.707107*const_cord_y;
}  
  
  
  
void calculos()
{
  float bar_v[3], bar_u[3];  /* baricentro do modelo e dados */
  float trans_mod[4][4], trans_dados[4][4],
        scal[4][4],  /* mat translacao e escala */
        aux1[4][4], aux2[4][4], aux1_inv[4][4], /* auxiliares */
        rot[4][4], mat_final[4][4]; /* rotacao e final */
  
  float var_v, var_u, alfa; /* variancia dos vert, dados e alfa */
  float aux4, aux5, aux3;
  float vet_aux1[3], vet_aux2[3];

  float d, dmin; /* para calculo de distancia */
  int indmin;  /* indice de menor distancia */

  float delta, discrepancia, valor;
  int ncasados;
  
  unsigned char teste;
  
  float normal_v1[3], normal_v2[3], normal_v3[3];/* normais dos vert.*/
  float normal_u1[3], normal_u2[3], normal_u3[3];/* normais dos dados*/
  
  int i, j, k,  /* indices utilizados p/ vertices do modelo */
      l, m, n;  /* indices utilizados p/ vertices dos dados */

  int t, h;
  
  for(i=0; i<nvertice_mod; i++)
    for(j=i+1; j<nvertice_mod; j++)
      for(k=j+1; k<nvertice_mod; k++)
	{
	  for(l=0; l<nvertice_dados; l++)
	    for(m=0; m<nvertice_dados; m++)
	      for(n=0; n<nvertice_dados; n++)
		if((l!=m)&&(l!=n)&&(m!=n))
		  {
		    /* ----- calculo do baricentro --- */
		    bar_v[0] = (mv[i][0] + mv[j][0] + mv[k][0])/3.0;
		    bar_v[1] = (mv[i][1] + mv[j][1] + mv[k][1])/3.0;
		    bar_v[2] = (mv[i][2] + mv[j][2] + mv[k][2])/3.0;
		    
		    bar_u[0] = (mu_3d[l][0] + mu_3d[m][0] + mu_3d[n][0])/3.0;
		    bar_u[1] = (mu_3d[l][1] + mu_3d[m][1] + mu_3d[n][1])/3.0;
		    bar_u[2] = (mu_3d[l][2] + mu_3d[m][2] + mu_3d[n][2])/3.0;

		    
		    
		    /* ------ Matriz de translacao --- */
		    trans_mod[0][0] = 1.0;
		    trans_mod[0][1] = - bar_v[0];
		    trans_mod[0][2] = - bar_v[1];
		    trans_mod[0][3] = - bar_v[2];
		    
		    trans_mod[1][0] = 0.0;
		    trans_mod[1][1] = 1.0;
		    trans_mod[1][2] = 0.0;
		    trans_mod[1][3] = 0.0;
		    
		    trans_mod[2][0] = 0.0;
		    trans_mod[2][1] = 0.0;
		    trans_mod[2][2] = 1.0;
		    trans_mod[2][3] = 0.0;
		    
		    trans_mod[3][0] = 0.0;
		    trans_mod[3][1] = 0.0;
		    trans_mod[3][2] = 0.0;
		    trans_mod[3][3] = 1.0;

		    trans_dados[0][0] = 1.0;
		    trans_dados[0][1] = bar_u[0] ;
		    trans_dados[0][2] = bar_u[1] ;
		    trans_dados[0][3] = bar_u[2] ;
		    
		    trans_dados[1][0] = 0.0;
		    trans_dados[1][1] = 1.0;
		    trans_dados[1][2] = 0.0;
		    trans_dados[1][3] = 0.0;
		    
		    trans_dados[2][0] = 0.0;
		    trans_dados[2][1] = 0.0;
		    trans_dados[2][2] = 1.0;
		    trans_dados[2][3] = 0.0;
		    
		    trans_dados[3][0] = 0.0;
		    trans_dados[3][1] = 0.0;
		    trans_dados[3][2] = 0.0;
		    trans_dados[3][3] = 1.0;

		    
		    
		    
	       /* ---------- Variancia total ---- */
		    aux4 = dist(mv[i],bar_v);
		    aux4 *= aux4;
		    aux5 = dist(mv[j],bar_v);
		    aux5 *= aux5;
		    aux3 = dist(mv[k],bar_v);
		    aux3 *= aux3;
		    var_v = aux4 + aux5 + aux3;
		    
		    aux4 = dist(mu_3d[l],bar_u);
		    aux4 *= aux4;
		    aux5 = dist(mu_3d[m],bar_u);
		    aux5 *= aux5;
		    aux3 = dist(mu_3d[n],bar_u);
		    aux3 *= aux3;
		    var_u = aux4 + aux5 + aux3;
		    
		    alfa = sqrt(var_u/var_v);

	       /* ---------- Matriz de escala --- */
		    scal[0][0] = 1.0;
		    scal[0][1] = 0.0;
		    scal[0][2] = 0.0;
		    scal[0][3] = 0.0;
		    scal[1][0] = 0.0;
		    scal[1][1] = alfa;
		    scal[1][2] = 0.0;
		    scal[1][3] = 0.0;
		    scal[2][0] = 0.0;
		    scal[2][1] = 0.0;
		    scal[2][2] = alfa;
		    scal[2][3] = 0.0;
		    scal[3][0] = 0.0;
		    scal[3][1] = 0.0;
		    scal[3][2] = 0.0;
		    scal[3][3] = alfa;
		    


	        /*----- --- calculo das normais --- */
		    vet_aux1[0] = mv[j][0] - mv[i][0];
		    vet_aux1[1] = mv[j][1] - mv[i][1];
		    vet_aux1[2] = mv[j][2] - mv[i][2];
		    
		    vet_aux2[0] = mv[k][0] - mv[i][0];
		    vet_aux2[1] = mv[k][1] - mv[i][1];
		    vet_aux2[2] = mv[k][2] - mv[i][2];
		    
		    crossp(vet_aux1, vet_aux2, normal_v1);

		    vet_aux1[0] = mu_3d[m][0] - mu_3d[l][0];
		    vet_aux1[1] = mu_3d[m][1] - mu_3d[l][1];
		    vet_aux1[2] = mu_3d[m][2] - mu_3d[l][2];

		    vet_aux2[0] = mu_3d[n][0] - mu_3d[l][0];
		    vet_aux2[1] = mu_3d[n][1] - mu_3d[l][1];
		    vet_aux2[2] = mu_3d[n][2] - mu_3d[l][2];
		    
		    crossp(vet_aux1, vet_aux2, normal_u1);

		/* ------- calculo de vetor ortogonal a normal_1 */
		    normal_v2[0] = mv[i][0] - bar_v[0];
		    normal_v2[1] = mv[i][1] - bar_v[1];
		    normal_v2[2] = mv[i][2] - bar_v[2];
		    
		    normalize(normal_v2);
		    
		    normal_u2[0] = mu_3d[l][0] - bar_u[0];
		    normal_u2[1] = mu_3d[l][1] - bar_u[1];
		    normal_u2[2] = mu_3d[l][2] - bar_u[2];
		    
		    normalize(normal_u2);
		    
	     /* ------ calculo do terceiro vetor ortogonal */
		    crossp(normal_v1, normal_v2, normal_v3);
		    crossp(normal_u1, normal_u2, normal_u3);


		/*----- matriz de rotacao ----- */

		    aux1[0][0] = 1.0;
		    aux1[0][1] = 0.0;
		    aux1[0][2] = 0.0;
		    aux1[0][3] = 0.0;
		    aux1[1][0] = 0.0;
		    aux1[1][1] = normal_v1[0];
		    aux1[1][2] = normal_v1[1];
		    aux1[1][3] = normal_v1[2];
		    aux1[2][0] = 0.0;
		    aux1[2][1] = normal_v2[0];
		    aux1[2][2] = normal_v2[1];
		    aux1[2][3] = normal_v2[2];
		    aux1[3][0] = 0.0;
		    aux1[3][1] = normal_v3[0];
		    aux1[3][2] = normal_v3[1];
		    aux1[3][3] = normal_v3[2];
		    
		    aux2[0][0] = 1.0;
		    aux2[0][1] = 0.0;
		    aux2[0][2] = 0.0;
		    aux2[0][3] = 0.0;
		    aux2[1][0] = 0.0;
		    aux2[1][1] = normal_u1[0];
		    aux2[1][2] = normal_u1[1];
		    aux2[1][3] = normal_u1[2];
		    aux2[2][0] = 0.0;
		    aux2[2][1] = normal_u2[0];
		    aux2[2][2] = normal_u2[1];
		    aux2[2][3] = normal_u2[2];
		    aux2[3][0] = 0.0;
		    aux2[3][1] = normal_u3[0];
		    aux2[3][2] = normal_u3[1];
		    aux2[3][3] = normal_u3[2];
		    
		    inversa(aux1, aux1_inv); /* calcula matriz
						inversa de aux1 */
		    
		    mat_prod(aux1_inv, aux2, rot); /* rot */



	       /* --------- matriz  final---------  */
		      
		    mat_prod(trans_mod, scal, aux1); /* aux1 = trans_mod*scal*/
	            mat_prod(aux1, rot, aux2);  /* aux2 = trans_mod*scal*rot */
                    mat_prod(aux2, trans_dados, mat_final);  
		      /* mat_final = trans_mod*scal*rot*trans_dados    */


	      /* ------ QUALIDADE DA MATRIZ FINAL ----*/
	      /* ----- vertices do modelo mapeados em mv_map_hom[4] -----*/

		    for(h=0; h<nvertice_mod; h++)
		      prod_row_mat(mv_homog[h], mat_final, mv_map_hom[h]);
		    
		    

		/*--- procura do mv_map_hom mais proximo de u ----*/
		
		    		
		    /* transforma vertices mapeados de 3d para isotropicos */
		    for(h=0; h<nvertice_mod; h++)
		      conv_3d_iso(mv_map_hom[h][1], mv_map_hom[h][2], mv_map_hom[h][3], &mv_map_iso[h][0], &mv_map_iso[h][1], &mv_map_iso[h][2]);

		    nvisiveis = 0;
		    ncasados = 0;
		    discrepancia =0.0;
		    for(h=0; h<nvertice_mod; h++)
		      {
			dmin = dist(mu_iso[0], mv_map_iso[h]);
			
			for(t=1; t<nvertice_dados; t++)
			  {
			    d = dist(mu_iso[t], mv_map_iso[h]);
			    if(d < dmin)
			      dmin = d;
			  }
			
			delta = dmin*dmin;
	/*		fprintf(stderr,"DELTA = %f\n",delta); */
			if(delta<deltamax)
			  {
			    ncasados++;
			    discrepancia += delta;
			  }
		      }

		  /* --- imprimir resultado -----*/

		    if(ncasados>=NUM_CASA)
		      {
			valor = sqrt(discrepancia/(ncasados + 1));
			fprintf(file_log," Ncasados = %d\n    Valor = %f\n",
				ncasados, valor);
			fprintf(file_log," Discrepancia = %f\n",discrepancia);
			fprintf(file_log," Vertices usados com modelo = [%d] [%d] [%d]\n",i,j,k);
			fprintf(file_log," Vertices usados com dados  = [%d] [%d] [%d]\n",l,m,n);
			
		/*	  fprintf(file_log,"  Matriz trans\n");
			  ver2(trans,4,4);
			  fprintf(file_log,"  Matriz Scal\n");
			  ver2(scal,4,4);
			  fprintf(file_log,"  Matriz rot\n");
			  ver2(rot,4,4);
			  fprintf(file_log,"  Matriz final\n");
			  ver2(mat_final,4,4);
			  fprintf(file_log,"  Matriz mu_3d dos dados\n");
			  ver1(mu_3d,7,3);
			  fprintf(file_log,"  Matriz mv dos vertices\n");
			  ver1(mv,8,3);   */
			  
			  fprintf(file_log,"  Matriz mv_map_hom[8][4] vertices mapeados\n");
			  ver1(mv_map_hom,nvertice_mod,4);

			  /* ESCRITA DOS PTOS MAPEADOS EM ARQUIVO
			     PARA PLOTAGEM EM CANVAS */

			fprintf(saida,"%9.4f %9.4f %9.4f   %3d %9.4f  %9.4f\n",mv_map_hom[0][1], mv_map_hom[0][2],mv_map_hom[0][3],ncasados,discrepancia, valor);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[1][1],
				mv_map_hom[1][2],mv_map_hom[1][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[0][1],
				mv_map_hom[0][2],mv_map_hom[0][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[3][1],
				mv_map_hom[3][2],mv_map_hom[3][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[0][1],
				mv_map_hom[0][2],mv_map_hom[0][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[6][1],
				mv_map_hom[6][2],mv_map_hom[6][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[1][1],
				mv_map_hom[1][2],mv_map_hom[1][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[2][1],
				mv_map_hom[2][2],mv_map_hom[2][3]);
			  
           /* */	fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[1][1],
				  mv_map_hom[1][2],mv_map_hom[1][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[7][1],
				  mv_map_hom[7][2],mv_map_hom[7][3]);
	      /* */	  
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[2][1],
				mv_map_hom[2][2],mv_map_hom[2][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[3][1],
				mv_map_hom[3][2],mv_map_hom[3][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[2][1],
				mv_map_hom[2][2],mv_map_hom[2][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[5][1],
				mv_map_hom[5][2],mv_map_hom[5][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[3][1],
				mv_map_hom[3][2],mv_map_hom[3][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[4][1],
				mv_map_hom[4][2],mv_map_hom[4][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[4][1],
				mv_map_hom[4][2],mv_map_hom[4][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[5][1],
				mv_map_hom[5][2],mv_map_hom[5][3]);
			
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[4][1],
				mv_map_hom[4][2],mv_map_hom[4][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[6][1],
				mv_map_hom[6][2],mv_map_hom[6][3]);
			
		/* */   fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[5][1],
				mv_map_hom[5][2],mv_map_hom[5][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[7][1],
				mv_map_hom[7][2],mv_map_hom[7][3]);
	     /*	*/	  
		/* */	fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[6][1],
				mv_map_hom[6][2],mv_map_hom[6][3]);
			fprintf(saida,"%9.4f %9.4f %9.4f\n",mv_map_hom[7][1],
				mv_map_hom[7][2],mv_map_hom[7][3]);
	     /*	*/		  
			fprintf(saida,"%f %f %f\n",-1.0,-1.0,-1.0);
			  
			  
			  
		      }
		    
		    
		  }
        }
}

void main(argc, argv)
int argc;
char *argv[];
{
  int i;
  char aux[255], carac;
  char parametros[50], saida_aux[50],saida_log[50], saida_mod[50],
       entrada_vert[50], modelo[50];
  
  double start, stop; /* para medida de tempo de execucao */
  FILE *arq_parametros, *file_modelo, *entrada;

  strcpy(parametros, argv[2]);
  strcat(parametros, ".parms");
  
  strcpy(saida_aux, argv[1]);
  strcat(saida_aux, "-");
  strcat(saida_aux, argv[2]);

  strcpy(saida_log, saida_aux);
  strcat(saida_log, "_mod.log");  

  strcpy(saida_mod, saida_aux);
  strcat(saida_mod, "_3d.mod");  
  
  strcpy(entrada_vert, saida_aux);
  strcat(entrada_vert, "_iso.vert");

  entrada = fopen(entrada_vert, "r");
  saida = fopen(saida_mod, "w");
  file_log = fopen(saida_log, "w");
  arq_parametros = fopen(parametros, "r");
  
  fgets(aux, 254, arq_parametros);  /* leitura do cabecalho   */
  fgets(aux, 254, arq_parametros);  /* leitura de duas linhas do arquivo de   */
  fgets(aux, 254, arq_parametros);  /* entrada que e usada apenas em tese-seg */
  fscanf(arq_parametros,"foco = %f\n",&foco);
  fscanf(arq_parametros,"baseline = %f\n",&baseline);
  
  const_cord_x = (SIZEX-1.0)/2.0;
  const_cord_y = (SIZEY-1.0)/2.0;
    
  fprintf(stderr,"\n\n   Casamento do modelo com dados obtidos.\n");

  fprintf(stderr,"\n\n       Dados.\n");
  fprintf(stderr,"\n Arquivo de vertices de entrada = %s",entrada_vert);
    
  fprintf(stderr,"\n Arquivo saida(v-mod-mapeados)  = %s",saida_mod);

  fprintf(file_log,"\n\n   Casamento do modelo com dados obtidos.\n");

  fprintf(file_log,"\n\n       Dados.\n");
  fprintf(file_log,"\n Arquivo de vertices de entrada = %s",entrada_vert);
    
  fprintf(file_log,"\n Arquivo saida(v-mod-mapeados)  = %s\n",saida_mod);


  fprintf(stderr,"\n Delta max   = ");
  scanf("%f",&deltamax);

  nvertice_dados = 0;

  while((fgets(aux, 254, entrada)) != NULL)
    nvertice_dados++;

  fprintf(stderr,"nvertice_dados = %d\n",nvertice_dados);

    
  fprintf(stderr,"\n Modelo abstrato = ");
  scanf("%s",modelo); 
  file_modelo = fopen(modelo, "r");

  nvertice_mod = 0;

  carac = getc(file_modelo);
  while(carac == '#')
    {
      fgets(aux, 254,file_modelo);
      carac = getc(file_modelo);
    }
  
  ungetc(carac, file_modelo);
  
  
  while(fgets(aux, 254, file_modelo) != NULL)
    nvertice_mod++;
	
  
  fprintf(stderr,"nvertice_mod = %d\n",nvertice_mod);
  

  mv = (float **)malloc(nvertice_mod * sizeof(float *));
  for(i=0; i<nvertice_mod; i++)
    mv[i] = (float *)malloc(3 * sizeof(float));

  mv_homog = (float **)malloc(nvertice_mod * sizeof(float *));
  for(i=0; i<nvertice_mod; i++)
    mv_homog[i] = (float *)malloc(4 * sizeof(float));

  mv_map_hom = (float **)malloc(nvertice_mod * sizeof(float *));
  for(i=0; i<nvertice_mod; i++)
    mv_map_hom[i] = (float *)malloc(4 * sizeof(float));

  mv_map_iso = (float **)malloc(nvertice_mod * sizeof(float *));
  for(i=0; i<nvertice_mod; i++)
    mv_map_iso[i] = (float *)malloc(3 * sizeof(float));


  mu_3d = (float **)malloc(nvertice_dados * sizeof(float *));
  for(i=0; i<nvertice_dados; i++)
    mu_3d[i] = (float *)malloc(3 * sizeof(float));


  mu_iso = (float **)malloc(nvertice_dados * sizeof(float *));
  for(i=0; i<nvertice_dados; i++)
    mu_iso[i] = (float *)malloc(3 * sizeof(float));

  rewind(entrada);

  for(i=0; i<nvertice_dados; i++)
    fscanf(entrada,"%f %f %f\n",&mu_iso[i][0],&mu_iso[i][1],&mu_iso[i][2]);

 /* transforma vertices isotropicos em 3d */
  for(i=0; i<nvertice_dados; i++)
    conv_iso_3d(mu_iso[i][0], mu_iso[i][1], mu_iso[i][2], &mu_3d[i][0], &mu_3d[i][1], &mu_3d[i][2]); 

    
/*  --- Leitura dos vertices do modelo ------------- */

  rewind(file_modelo);
  
  carac = getc(file_modelo);
  while(carac == '#')
    {
      fgets(aux, 254,file_modelo);
      carac = getc(file_modelo);
    }
  
  ungetc(carac, file_modelo);

  for(i=0; i<nvertice_mod; i++)
    fscanf(file_modelo,"%f %f %f\n",&mv[i][0],&mv[i][1],&mv[i][2]);

/* vertices do modelo em coordenadas homogeneas :    */    

  for(i=0; i<nvertice_mod; i++)
    {
      mv_homog[i][0] = 1.0;
      mv_homog[i][1] = mv[i][0];
      mv_homog[i][2] = mv[i][1];
      mv_homog[i][3] = mv[i][2];
    }

    fprintf(file_log,"  Matriz mu_iso dos dados\n");
    ver1(mu_iso,nvertice_dados,3);

    fprintf(file_log,"  Matriz mu_3d dos dados\n");    
    ver1(mu_3d,nvertice_dados,3);
    
    fprintf(file_log,"  Matriz mv dos vertices\n");
    ver1(mv,nvertice_mod,3);
    
    fprintf(file_log,"  Matriz mv dos vertices homogeneos\n");
    ver1(mv_homog,nvertice_mod,4);
 
    calculos();
    
}