dist2isosurf_.c Source File

00001 
00015 #include "interpolation.h"
00016 #include "utils.h"
00017 #include <inrimage/image.h>
00018 
00019 extern void imerror ( int code, char * format, ...);
00020 float mediane(char * famille,int j,int buff_length,int dim_x,Fort_int nb_octets,float next);
00021 
00033 void dist2isosurf(char *distances,char * famille,float marge,int nb_fam,int time_flag)
00034 {
00035   struct nf_fmt gfmt_in,gfmt_out;
00036 
00037   struct image * img_distances;
00038   struct image * img_famille;
00039 
00040   //  const char* format_table[3]={"REELLE","FIXE","PACKEE"};
00041   char ** buff_famille=NULL;
00042 
00043   clock_t t0=0,t1=0;
00044 
00045   int i,k,z_in,z_out,nv_zero;
00046 
00047   /* chronometrage si on le demande */
00048   if (time_flag)
00049     t0 = clock();
00050 
00051   /* init */
00052 
00053   for (i=0;i<9;i++)
00054     gfmt_in.lfmt[i]=0;
00055   for (i=0;i<3;i++)
00056     gfmt_in.offsets[i]=0;
00057   gfmt_in.maille=0;
00058   gfmt_in.bias=0;
00059   gfmt_in.scale=1;
00060 
00061   /* ouverture des images a convertir, et recuperation de ses parametres */
00062   img_distances=image_(distances,"e","",&gfmt_in);
00063 
00064   /* verification de la valeur de Z */
00065   if ((gfmt_in.lfmt[I_NDIMZ]<2))
00066     imerror(6,"dist2isosurf: Erreur: image en 3D requise pour %s, valeur lue: %d\n",distances,gfmt_in.lfmt[I_NDIMZ]);
00067 
00068   z_in = gfmt_in.lfmt[I_NDIMZ];
00069 
00070   buff_famille=separation(img_distances,gfmt_in,&gfmt_out,marge,nb_fam,&nv_zero,max(time_flag-1,0));
00071 
00072   z_out = gfmt_out.lfmt[I_NDIMZ];
00073   /* creation de l'image de sortie */
00074   img_famille=image_(famille,"c","a",&gfmt_out);
00075 
00076   /* ecriture de la famille de surfaces */
00077 
00078   switch(gfmt_out.lfmt[I_BSIZE])
00079     {
00080     case 1:
00081       for (k=0;k<z_out;k++)
00082          ecr_(&img_famille,&(gfmt_in.lfmt[I_NDIMY]),buff_famille[k]);
00083       break;
00084     case 4:
00085       for (k=0;k<z_out;k++)
00086          ecrflt_(&img_famille,&(gfmt_in.lfmt[I_NDIMY]),buff_famille[k]);
00087       break;
00088     }
00089   /* liberation des variables temporaires et fermeture des images */
00090 
00091   for (k=0;k<z_out;k++)
00092     i_Free(&(buff_famille[k]));
00093   free(buff_famille);
00094   buff_famille=NULL;
00095 
00096   fermnf_(&img_distances);
00097   fermnf_(&img_famille);
00098 
00099   if (time_flag)
00100     {
00101       t1=clock()-t0;
00102       printf("dist2isosurf: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
00103     }
00104 }
00105 
00121 char ** separation(struct image * img_distances,struct nf_fmt gfmt_in,struct nf_fmt * gfmt_out,float marge,int nb_fam,int * nv_zero,int time_flag)
00122 {
00123   char * buff_distances=NULL;
00124   char ** famille=NULL;
00125 
00126   float minDist=0,maxDist=0;
00127   float epsilon;
00128 
00129   int i,j,k;
00130   int z_in,z_out,dim_x;
00131   int buff_size,buff_length=0;
00132   int num_plan;
00133 
00134   clock_t t0=0,t1=0;
00135   
00136   /* chronometrage si on le demande */
00137   if (time_flag)
00138     t0 = clock();
00139 
00140   z_in = gfmt_in.lfmt[I_NDIMZ];
00141 
00142   /* initialisation de gfmt_out */
00143   for (i=0;i<9;i++)
00144     gfmt_out->lfmt[i]=gfmt_in.lfmt[i];
00145   for (i=0;i<3;i++)
00146     gfmt_out->offsets[i]=gfmt_in.offsets[i];
00147   gfmt_out->maille=gfmt_in.maille;
00148   gfmt_out->bias=gfmt_in.bias;
00149   gfmt_out->scale=gfmt_in.scale;
00150 
00151   /* calculer min et max de la distance */
00152 
00153   extrema3D(img_distances,gfmt_in.lfmt,&minDist,&maxDist);
00154 
00155   /* en deduire epsilon */
00156 
00157   if (nb_fam!=0)
00158     {
00159       /* on a un nombre de familles imposees, on en deduit epsilon */
00160       gfmt_out->lfmt[I_NDIMZ]=nb_fam;
00161     }
00162   else if ((maxDist-minDist)<marge) /* l'epsilon demande est trop gros */
00163     {
00164       gfmt_out->lfmt[I_NDIMZ]=1;
00165     }
00166   else /* l'epsilon permet de calculer le nombre de surfaces de la famille, on le calcule */
00167     {
00168       gfmt_out->lfmt[I_NDIMZ]=(int)((maxDist-minDist)/marge);
00169     }
00170 
00171   gfmt_out->lfmt[I_DIMY]=gfmt_out->lfmt[I_NDIMY]*gfmt_out->lfmt[I_NDIMZ];
00172 
00173   z_out = gfmt_out->lfmt[I_NDIMZ];
00174   epsilon=(maxDist-minDist)/(float)(z_out);
00175 
00176   famille=(char **)malloc (z_out*sizeof(char*));
00177 
00178   /* pour chaque plan remplir de 0 le resultat */
00179   buff_size =  buffer_size(gfmt_out->lfmt,gfmt_out->lfmt[I_NDIMY]);
00180 
00181   switch(gfmt_in.lfmt[I_BSIZE])
00182     {
00183     case 1: /* image de char */
00184       buff_length = buff_size/sizeof(char);
00185       
00186       /* buff_famille represente une famille de surfaces, chaque element de la famille correspondant a un plan */
00187       for(k=0;k<z_out;k++)
00188          {
00189            famille[k]=i_malloc(buff_size);
00190            for (j=0;j<buff_length;j++)
00191              famille[k][j]=0;
00192          }
00193       break;
00194     case 4: /* image de float */
00195       buff_length = buff_size/sizeof(float);
00196       
00197       /* buff_famille represente une famille de surfaces, chaque element de la famille correspondant a un plan */      
00198       for(k=0;k<z_out;k++)
00199          {
00200            famille[k]=i_malloc(buff_size);
00201            for (j=0;j<buff_length;j++)
00202              ((float *)(famille[k]))[j]=0;
00203          }
00204       break;
00205     default:
00206       imerror(6,"separation: format non reconnu, codage ni sur 1 ni sur 4 octets\n");
00207       break;
00208     }
00209 
00210   buff_distances = i_malloc(buff_size);
00211 
00212   *nv_zero=(int)min(z_out-1,(-minDist)/epsilon);
00213 
00214   /* charger chaque plan et placer la valeur dans le bon plan */
00215 
00216   switch(gfmt_in.lfmt[I_BSIZE])
00217     {
00218     case 1: /* image de char */
00219       for (k=0;k<z_in;k++)
00220          {
00221            lect_(&img_distances,&(gfmt_out->lfmt[I_NDIMY]),buff_distances);
00222            for (j=0;j<buff_length;j++)
00223              {
00224                num_plan = (int)min(z_out-1,(((float)(buff_distances[j])-minDist)/epsilon));
00225                famille[num_plan][j]=(char)max((float)(famille[num_plan][j]),255*(float)k/(float)z_in);
00226              }
00227          }
00228       break;
00229     case 4: /* image de float */
00230       for (k=0;k<z_in;k++)
00231          {
00232            lect_(&img_distances,&(gfmt_out->lfmt[I_NDIMY]),buff_distances);
00233            for (j=0;j<buff_length;j++)
00234              {
00235                num_plan = (int)min(z_out-1,((((float *)buff_distances)[j]-minDist)/epsilon));
00236                ((float *)(famille[num_plan]))[j]=max(((float *)(famille[num_plan]))[j],(float)k/(float)z_in);
00237              }
00238          }
00239       break;
00240     }
00241   /* lisser en partant de l'isosurface de niveau zero */
00242   
00243   dim_x=gfmt_out->lfmt[I_NDIMX];
00244 
00245   switch(gfmt_in.lfmt[I_BSIZE])
00246     {
00247     case 1: /* image de caracteres */
00248       for (k=*nv_zero;k<z_out;k++)
00249          /* on lisse en se referant aux valeurs autour et a la valeur au meme point de la surface precedente */
00250          for (j=0;j<buff_length;j++)
00251            if (famille[k][j]==0)
00252              famille[k][j]=(char)mediane(famille[k],j,buff_length,dim_x,gfmt_out->lfmt[I_BSIZE],(float)(famille[(int)max(k-1,0)][j]));       
00253       
00254       for (k=*nv_zero-1;k>=0;k--)
00255          /* on lisse en se referant aux valeurs autour et a la valeur au meme point de la surface suivante */
00256          for (j=0;j<buff_length;j++)
00257            if (famille[k][j]==0)
00258              famille[k][j]=(char)mediane(famille[k],j,buff_length,dim_x,gfmt_out->lfmt[I_BSIZE],(float)(famille[(int)min(k+1,z_out)][j]));
00259       break;
00260     case 4: /* image de float */
00261       for (k=*nv_zero;k<z_out;k++)
00262          /* on lisse en se referant aux valeurs autour et a la valeur au meme point de la surface precedente */
00263          for (j=0;j<buff_length;j++)
00264            if (((float *)(famille[k]))[j]==0)
00265              ((float *)(famille[k]))[j]=mediane(famille[k],j,buff_length,dim_x,gfmt_out->lfmt[I_BSIZE],((float *)(famille[(int)max(k-1,0)]))[j]);
00266       
00267       for (k=*nv_zero-1;k>=0;k--)
00268          /* on lisse en se referant aux valeurs autour et a la valeur au meme point de la surface suivante */
00269          for (j=0;j<buff_length;j++)
00270            if (((float *)(famille[k]))[j]==0)
00271              ((float *)(famille[k]))[j]=mediane(famille[k],j,buff_length,dim_x,gfmt_out->lfmt[I_BSIZE],((float *)(famille[(int)min(k+1,z_out)]))[j]);
00272       break;
00273     }
00274   i_Free(&buff_distances);
00275   if (time_flag)
00276     {
00277       t1=clock()-t0;
00278       printf("separation: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
00279     }
00280 
00281   return famille;
00282 }
00283 
00284 float moyenne(float a,float b,float c,float d,float e,float f,float g,float h,float i);
00285 
00286 float mediane(char * famille,int j,int buff_length,int dim_x,Fort_int nb_octets,float next)
00287 {
00288   /* premiere ligne */
00289   
00290   if (j<dim_x)
00291     {
00292       /* premiere case premiere ligne */
00293       if (j == 0)
00294          switch(nb_octets)
00295            {
00296            case 1:
00297              return moyenne(next,famille[1],famille[dim_x],famille[dim_x+1],0,0,0,0,0);
00298              /*       return ((next + famille[1] +
00299                        famille[dim_x] + famille[dim_x+1])
00300                        /4);*/ break;
00301            case 4:
00302              return moyenne(next,((float *)famille)[1],((float *)famille)[dim_x],((float *)famille)[dim_x+1],0,0,0,0,0);          
00303              /*       return ((next + ((float *)famille)[1] +
00304                        ((float *)famille)[dim_x] + ((float *)famille)[dim_x+1])
00305                        /4);*/ break;
00306            }
00307       /* derniere case premiere ligne */
00308       if (j + 1 == dim_x ) 
00309          switch(nb_octets)
00310            {
00311            case 1:
00312              return moyenne(next,famille[j-1],famille[j+dim_x],famille[j+dim_x-1],0,0,0,0,0);
00313              /*       return ((next + famille[j-1] +
00314                        famille[j+dim_x] + famille[j+dim_x-1])
00315                        /4);*/ break;
00316            case 4:
00317              return moyenne(next,((float *)famille)[j-1],((float *)famille)[j+dim_x],((float *)famille)[j+dim_x-1],0,0,0,0,0);
00318              /*       return ((next + ((float *)famille)[j-1] +
00319                        ((float *)famille)[j+dim_x] + ((float *)famille)[j+dim_x-1])
00320                        /4);*/ break;
00321            }
00322       /* premiere ligne */
00323       switch(nb_octets)
00324          {
00325          case 1:
00326            return moyenne(next,famille[j-1],famille[j+1],famille[j+dim_x-1],famille[j+dim_x],famille[j+dim_x+1],0,0,0);
00327            /*       return ((next + famille[j-1] + famille[j+1] +
00328                      famille[j+dim_x-1] + famille[j+dim_x] + famille[j+dim_x+1])
00329                      /6);*/break;
00330          case 4:
00331            return moyenne(next,((float *)famille)[j-1],((float *)famille)[j+1],((float *)famille)[j+dim_x-1],((float *)famille)[j+dim_x],((float *)famille)[j+dim_x+1],0,0,0);
00332            /*       return ((next + ((float *)famille)[j-1] + ((float *)famille)[j+1] +
00333                      ((float *)famille)[j+dim_x-1] + ((float *)famille)[j+dim_x] + ((float *)famille)[j+dim_x+1])
00334                      /6);*/break;
00335          }
00336     }
00337 
00338   /* derniere ligne */
00339 
00340   if (j>buff_length-dim_x)
00341     {
00342       /* premiere case derniere ligne */
00343       if (j%dim_x == 0)
00344          switch(nb_octets)
00345            {
00346            case 1:
00347              return moyenne(next,famille[j+1],famille[j-dim_x],famille[j-dim_x+1],0,0,0,0,0);
00348              /*       return ((next + famille[j+1] +
00349                        famille[j-dim_x] + famille[j-dim_x+1])
00350                        /4);*/ break;
00351            case 4:
00352              return moyenne(next,((float *)famille)[j+1],((float *)famille)[j-dim_x],((float *)famille)[j-dim_x+1],0,0,0,0,0);
00353              /*       return ((next + ((float *)famille)[j+1] +
00354                        ((float *)famille)[j-dim_x] + ((float *)famille)[j-dim_x+1])
00355                        /4);*/ break;
00356            }
00357  
00358       /* derniere case derniere ligne */
00359       if ((j + 1) % dim_x == 0 )
00360          switch(nb_octets)
00361            {
00362            case 1:
00363              return moyenne(next,famille[j-1],famille[j-1-dim_x],famille[j-dim_x],0,0,0,0,0);
00364              /*       return ((next + famille[j-1] +
00365                        famille[j-1-dim_x] + famille[j-dim_x])
00366                        /4);*/ break;
00367            case 4:
00368              return moyenne(next,((float *)famille)[j-1],((float *)famille)[j-1-dim_x],((float *)famille)[j-dim_x],0,0,0,0,0);
00369              /*       return ((next + ((float *)famille)[j-1] +
00370                        ((float *)famille)[j-1-dim_x] + ((float *)famille)[j-dim_x])
00371                        /4);*/ break;
00372            }
00373       /* derniere ligne */
00374       switch(nb_octets)
00375          {
00376          case 1:
00377            return moyenne(next,famille[j-1],famille[j+1],famille[j-1-dim_x],famille[j-dim_x],famille[j+1-dim_x],0,0,0);
00378            /*       return ((next + famille[j-1] + famille[j+1] +
00379                      famille[j-1-dim_x] + famille[j-dim_x] + famille[j+1-dim_x])
00380                      /6);*/ break;
00381          case 4:
00382            return moyenne(next,((float *)famille)[j-1],((float *)famille)[j+1],((float *)famille)[j-1-dim_x],((float *)famille)[j-dim_x],((float *)famille)[j+1-dim_x],0,0,0);
00383            /*       return ((next + ((float *)famille)[j-1] + ((float *)famille)[j+1] +
00384                      ((float *)famille)[j-1-dim_x] + ((float *)famille)[j-dim_x] + ((float *)famille)[j+1-dim_x])
00385                      /6);*/ break;
00386          }
00387     }
00388 
00389   /* ligne quelconque */
00390 
00391   /* premiere case */
00392   if (j%dim_x == 0)
00393     switch(nb_octets)
00394       {
00395       case 1: 
00396          return moyenne(famille[j-dim_x],famille[j-dim_x+1],next,famille[j+1],famille[j+dim_x],famille[j+dim_x+1],0,0,0);
00397          /*       return ((famille[j-dim_x] + famille[j-dim_x+1] +
00398                    next + famille[j+1] +
00399                    famille[j+dim_x] + famille[j+dim_x+1])
00400                    /6);*/ break;
00401       case 4:
00402          return moyenne(((float *)famille)[j-dim_x],((float *)famille)[j-dim_x+1],next,((float *)famille)[j+1],((float *)famille)[j+dim_x],((float *)famille)[j+dim_x+1],0,0,0);
00403          /*       return ((((float *)famille)[j-dim_x] + ((float *)famille)[j-dim_x+1] +
00404                    next + ((float *)famille)[j+1] +
00405                    ((float *)famille)[j+dim_x] + ((float *)famille)[j+dim_x+1])
00406                    /6);*/ break;
00407       }
00408   /* derniere case */
00409   if ((j+1)%dim_x == 0)
00410     switch(nb_octets)
00411       {
00412       case 1: 
00413          return moyenne(famille[j-dim_x-1],famille[j-dim_x],famille[j-1],next,famille[j-1+dim_x],famille[j+dim_x],0,0,0);
00414          /*       return ((famille[j-dim_x-1] + famille[j-dim_x] +
00415                   famille[j-1] + next +
00416                   famille[j-1+dim_x] + famille[j+dim_x])
00417                   /6);*/ break;
00418       case 4: 
00419          return moyenne(((float *)famille)[j-dim_x-1],((float *)famille)[j-dim_x],((float *)famille)[j-1],next,((float *)famille)[j-1+dim_x],((float *)famille)[j+dim_x],0,0,0);
00420          /*       return ((((float *)famille)[j-dim_x-1] + ((float *)famille)[j-dim_x] +
00421                    ((float *)famille)[j-1] + next +
00422                    ((float *)famille)[j-1+dim_x] + ((float *)famille)[j+dim_x])
00423                    /6);*/ break;
00424       }
00425   switch(nb_octets)
00426     {
00427     case 1: 
00428       return moyenne(famille[j-dim_x-1],famille[j-dim_x],famille[j-dim_x+1],famille[j-1],next,famille[j+1],famille[j+dim_x-1],famille[j+dim_x],famille[j+dim_x+1]);
00429       /*      return ((famille[j-dim_x-1] + famille[j-dim_x] + famille[j-dim_x+1] +
00430                 famille[j-1] + next + famille[j+1] +
00431                 famille[j+dim_x-1] + famille[j+dim_x] + famille[j+dim_x+1])
00432                 /9);*/ break;
00433     case 4:
00434       return moyenne(((float *)famille)[j-dim_x-1],((float *)famille)[j-dim_x],((float *)famille)[j-dim_x+1],((float *)famille)[j-1],next,((float *)famille)[j+1],((float *)famille)[j+dim_x-1],((float *)famille)[j+dim_x],((float *)famille)[j+dim_x+1]);
00435       /*      return ((((float *)famille)[j-dim_x-1] + ((float *)famille)[j-dim_x] + ((float *)famille)[j-dim_x+1] + 
00436                 ((float *)famille)[j-1] + next + ((float *)famille)[j+1] + 
00437                 ((float *)famille)[j+dim_x-1] + ((float *)famille)[j+dim_x] + ((float *)famille)[j+dim_x+1])
00438                 /9);*/ break;
00439     }
00440   return 0;
00441 }
00442 
00443 float moyenne(float a,float b,float c,float d,float e,float f,float g,float h,float i)
00444 {
00445   int nb_non_nuls=0;
00446   if (a!=0)
00447     nb_non_nuls++;
00448   if (b!=0)
00449     nb_non_nuls++;
00450   if (c!=0)
00451     nb_non_nuls++;
00452   if (d!=0)
00453     nb_non_nuls++;
00454   if (e!=0)
00455     nb_non_nuls++;
00456   if (f!=0)
00457     nb_non_nuls++;
00458   if (g!=0)
00459     nb_non_nuls++;
00460   if (h!=0)
00461     nb_non_nuls++;
00462   if (i!=0)
00463     nb_non_nuls++;
00464   if (nb_non_nuls)
00465     return ((a+b+c+d+e+f+g+h+i)/nb_non_nuls);
00466   else
00467     return (0);
00468 }