---

G3DC_.c

Aller à la documentation de ce fichier.

#include "interpolation.h"
#include <unistd.h>

extern void imerror(int code, char * format, ...);

void calcul_parametres_differentiels(inricache *W,inricache *A,inricache *B,inricache *C,inricache *D,inricache *E,inricache *F,Fort_int lfmt[9],int time_flag);
void calcul_evolution_surface(inricache *imgPhi,inricache *imgForces,inricache *imgGeod,inricache *imgA,inricache *imgB,inricache *imgC,inricache *imgD,inricache *imgE,inricache *imgF,int nb_iter,float ptTau,Fort_int lfmt[9],int time_flag);
void calcul_force(inricache *imgPhi,inricache *imgForces,inricache *imgA,inricache *imgB,inricache *imgC,inricache *imgD,inricache *imgE,inricache *imgF,Fort_int lfmt[9],int time_flag);
void update_solution(inricache *imgPhi,inricache *imgForces,inricache *imgGeod,int nb_iter,float ptTau,Fort_int lfmt[9],int time_flag);

void g3Dc(char * distSrc,char * tempoW,char * distGeod,int nb_iter,float ptTau,int taille_cache,int time_flag)
{
  inricache *imgSrc,*imgGeod,*imgW;
  
  /* images temporaires */
  inricache *imgPhi,*imgForces;
  /* parametres differentiels */
  inricache *imgA,*imgB,*imgC,*imgD,*imgE,*imgF;

  struct nf_fmt gfmt;

  clock_t t0=0,t1=0;

  /* chronometrage si on le demande */
  if (time_flag)
    t0 = clock();
  t1=t0;

  fprintf(stderr,"creation des images temporaires ...");
  
  /* chargement / creation des images necessaires au calcul */

  /* le format est le meme pour toutes les images du processus */
  imgW=cree_inricache(tempoW,"e","",&gfmt,taille_cache);
  imgSrc=cree_inricache(distSrc,"e","a",&gfmt,taille_cache);

  imgGeod=cree_inricache(distGeod,"c","a",&gfmt,taille_cache);

  /* images temporaires pour le calcul de phi et F */
  imgPhi=cree_inricache("phi_cache","c","a",&gfmt,taille_cache);
  imgForces=cree_inricache("tempoF_cache","c","a",&gfmt,taille_cache);

  /* images temporaires correspondant aux parametres differentiels de W */

  imgA=cree_inricache("A_cache","c","a",&gfmt,taille_cache);
  imgB=cree_inricache("B_cache","c","a",&gfmt,taille_cache);
  imgC=cree_inricache("C_cache","c","a",&gfmt,taille_cache);
  imgD=cree_inricache("D_cache","c","a",&gfmt,taille_cache);
  imgE=cree_inricache("E_cache","c","a",&gfmt,taille_cache);
  imgF=cree_inricache("F_cache","c","a",&gfmt,taille_cache);

  if (max(time_flag-1,0))
    {
      t1=clock()-t1;
      fprintf(stderr,"(%f secondes)\n",t1/(float)CLOCKS_PER_SEC);
    }

  fprintf(stderr,"initialisation des buffers ...");

  /* on initialise le buffer de phi et de Geod avec l'image distSrc */

  copie_inricache(imgSrc,imgPhi);
  copie_inricache(imgSrc,imgGeod);

  /* on n'a plus besoin d'imgSrc */
  libere_inricache(imgSrc);

  if (max(time_flag-1,0))
    {
      t1=clock()-t1;
      fprintf(stderr,"(%f secondes)\n",t1/(float)CLOCKS_PER_SEC);
    }

  /* calcul de A,B,C,D,E,F,p,q,r */
  fprintf(stderr,"calcul des parametres differentiels de W ...");

  calcul_parametres_differentiels(imgW,
                                      imgA,imgB,imgC,imgD,imgE,imgF,
                                      gfmt.lfmt,(int)max(time_flag-2,0));

  /* on n'a plus besoin de imgW */
  libere_inricache(imgW);

  if (max(time_flag-1,0))
    {
      t1=clock()-t1;
      fprintf(stderr,"(%f secondes)\n",t1/(float)CLOCKS_PER_SEC);
    }

  /* calcul de l'evolution de la surface */
  fprintf(stderr,"calcul de l'evolution de la surface ...");


  /* on n'a plus besoin d'utiliser W, les coefficients _A,_B,_C,_D,_E et _F suffisent pour le calcul des forces*/
  calcul_evolution_surface(imgPhi,imgForces,imgGeod,
                              imgA,imgB,imgC,imgD,imgE,imgF,
                              nb_iter,ptTau,gfmt.lfmt,(int)max(time_flag-2,0));

  if (max(time_flag-1,0)==1)
    {
      t1=clock()-t1;
      fprintf(stderr,"(%f secondes)\n",t1/(float)CLOCKS_PER_SEC);
    }

  /* fermeture et sauvegarde des images */
  
  libere_inricache(imgGeod);
  libere_inricache(imgA);
  libere_inricache(imgB);
  libere_inricache(imgC);
  libere_inricache(imgD);
  libere_inricache(imgE);
  libere_inricache(imgF);
  libere_inricache(imgForces);

  /* suppression Physique des images temporaires * /

  unlink("A_cache");
  unlink("B_cache");
  unlink("C_cache");
  unlink("D_cache");
  unlink("E_cache");
  unlink("F_cache");
  */
  if (time_flag)
    {
      t1 = clock()-t0;
      printf("g3Dc: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
    }
}

void calcul_parametres_differentiels(inricache *W,inricache *A,inricache *B,inricache *C,inricache *D,inricache *E,inricache *F,Fort_int lfmt[9],int time_flag)
{
  register int i,j,k;
  clock_t t0=0,t1=0;
  float tmp_p,tmp_q,tmp_r,q2,p2,r2,denominateur,wzm,wzp,wym,wyp,wxm,wxp;

  /* chronometrage si on le demande */
  if (time_flag)
    t0 = clock();

  //  fprintf(stderr,"max=(%d,%d,%d)",lfmt[I_NDIMX],lfmt[I_NDIMY],lfmt[I_NDIMZ]);
  for (k=1;k<lfmt[I_NDIMZ]-1;k++)
    for (j=1;j<lfmt[I_NDIMY]-1;j++)
      for (i=1;i<lfmt[I_NDIMX]-1;i++)
         {
           wzm = getValue(W,i,j,k-1);
           wym = getValue(W,i,j-1,k);
           wxm = getValue(W,i-1,j,k);
           wxp = getValue(W,i+1,j,k);
           wyp = getValue(W,i,j+1,k);
           wzp = getValue(W,i,j,k+1);

           /* W en lecture , A B C D E F en ecriture */
           /* derivee partielle en x en (i,j,k) */
           tmp_p=(wxp - wxm)/2.0;

           /* derivee partielle en y en (i,j,k) */
           tmp_q=(wyp - wym)/2.0;
           
           /* derivee partielle en z en (i,j,k) */
           tmp_r=(wzp - wzm)/2.0;
           
           p2=tmp_p*tmp_p;
           q2=tmp_q*tmp_q;
           r2=tmp_r*tmp_r;
           denominateur=1+p2+q2+r2;
           
           setValue(A,i,j,k,(1+q2+r2)/denominateur);
           setValue(B,i,j,k,(1+p2+r2)/denominateur);
           setValue(C,i,j,k,(1+p2+q2)/denominateur);
           
           setValue(D,i,j,k,(2*tmp_p*tmp_q)/denominateur);
           setValue(E,i,j,k,(2*tmp_p*tmp_r)/denominateur);
           setValue(F,i,j,k,(2*tmp_q*tmp_r)/denominateur);

         }

  if (time_flag)
    {
      t1 = clock()-t0;
      printf("calcul_parametres_differentiels: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
    }
}

void calcul_evolution_surface(inricache *imgPhi,inricache *imgForces,inricache *imgGeod,inricache *imgA,inricache *imgB,inricache *imgC,inricache *imgD,inricache *imgE,inricache *imgF,int nb_iter,float ptTau,Fort_int lfmt[9],int time_flag)
{
  int i;
  clock_t init=0,prev=0,curr=0;

  /* chronometrage si on le demande. ici on ne chronometrera que le temps d'execution des iterations */
  if (time_flag)
    init = clock();
  curr=init;
  prev=init;
  
  for (i=0;i<nb_iter;i++)
    {
      calcul_force(imgPhi,imgForces,imgA,imgB,imgC,imgD,imgE,imgF,lfmt,(int)max(time_flag-1,0));
      update_solution(imgPhi,imgForces,imgGeod,ptTau,nb_iter,lfmt,(int)max(time_flag-1,0));

      if (time_flag)
         {
           curr=clock();
           fprintf(stderr,"calcul_evolution_surface: iteration %d,temps effectue: %f secondes,moyenne: %f secondes\n",i,(curr-prev)/(float)(CLOCKS_PER_SEC),((curr-init)/(float)(CLOCKS_PER_SEC)/(i+1)));
           prev=curr;
         }
    }
}

void calcul_force(inricache *imgPhi,inricache *imgForces,inricache *imgA,inricache *imgB,inricache *imgC,inricache *imgD,inricache *imgE,inricache *imgF,Fort_int lfmt[9],int time_flag)
{
  /* parametres differentiels */
  float a=0,b=0,c=0,d=0,e=0,f=0;
  /* derivees partielles en x,y,z de Phi dans les sens positifs et negatifs */
  float dxp=0,dxm=0,dyp=0,dym=0,dzp=0,dzm=0;
  /* derivees partielles en x,y,z de Phi par la methode du minmod */
  float dx=0,dy=0,dz=0;
  /* carre des derivees partielles en x,y,z */
  float dx2=0,dy2=0,dz2=0;
  float currPhi=0;

  register int i,j,k,i2,j2;

  clock_t t0=0,t1=0;

  /* chronometrage si on le demande */
  if (time_flag)
    t0=clock();

  /* todo: gestion des effets de bords, pixels extremaux */
  for (k=1;k<lfmt[I_NDIMZ]-1;k++)
    {
      for (j=1;j<lfmt[I_NDIMY]-1;j++)
         {
           for (i=1;i<lfmt[I_NDIMX]-1;i++)
             {

               a=getValue(imgA,i,j,k);
               b=getValue(imgB,i,j,k);
               c=getValue(imgC,i,j,k);
               d=getValue(imgD,i,j,k);
               e=getValue(imgE,i,j,k);
               f=getValue(imgF,i,j,k);

               currPhi = getValue(imgPhi,i,j,k);

               dxp = getValue(imgPhi,i+1,j,k) - currPhi;
               dxm = currPhi - getValue(imgPhi,i-1,j,k);

               dyp = getValue(imgPhi,i,j+1,k) - currPhi;
               dym = currPhi - getValue(imgPhi,i,j-1,k);

               dzp = getValue(imgPhi,i,j,k+1) - currPhi;
               dzm = currPhi - getValue(imgPhi,i,j,k-1);

               dx=minmod(dxp,dxm);
               dy=minmod(dyp,dym);
               dz=minmod(dzp,dzm);
               
               dx2=max3(dxp,-dxm,0);
               dx2*=dx2;
               dy2=max3(dyp,-dym,0);
               dy2*=dy2;
               dz2=max3(dzp,-dzm,0);           
               dz2*=dz2;

               setValue(imgForces,i,j,k,sqrt ((double)max(a*dx2+b*dy2+c*dz2-d*dx*dy-e*dx*dz-f*dy*dz,0)));
             }

           /* fin de la boucle en x prolongation des bords */
           for (i2=0;i2<lfmt[I_NDIMX];i2++)
             {
               setValue(imgForces,i2,0,k,getValue(imgForces,i2,1,k));
               setValue(imgForces,i2,lfmt[I_NDIMY]-1,k,getValue(imgForces,i2,lfmt[I_NDIMY]-2,k));
             }
           /* fin de la boucle en y prolongation des bords */
           for (j2=0;j2<lfmt[I_NDIMY];j2++)
             {
               setValue(imgForces,0,j2,k,getValue(imgForces,1,j2,k));
               setValue(imgForces,lfmt[I_NDIMX]-1,j2,k,getValue(imgForces,lfmt[I_NDIMX]-2,j2,k));
             }
         }
    }
  /*  for (j=0;j<lfmt[I_NDIMY];j++)
    for (i=0;i<lfmt[I_NDIMX];i++)
      {
         ((float *)(bufferF[0]))[indice(i,j,lfmt[I_NDIMX])] = ((float *)(bufferF[1]))[indice(i,j,lfmt[I_NDIMX])];
         ((float *)(bufferF[lfmt[I_NDIMZ]-1]))[indice(i,j,lfmt[I_NDIMX])] = ((float *)(bufferF[lfmt[I_NDIMZ]-2]))[indice(i,j,lfmt[I_NDIMX])];
      }
  */

  if (time_flag)
    {
      t1 = clock()-t0;
      printf("calcul_forces: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
    }

}

void update_solution(inricache *imgPhi,inricache *imgF,inricache *imgGeod,int nb_iter,float ptTau,Fort_int lfmt[9],int time_flag)
{
  register int i,j,k;
  float nvPhi,oldPhi;

  clock_t t0=0,t1=0;
  /* chronometrage si on le demande */
  if (time_flag)
    t0=clock();

  for (k=0;k<lfmt[I_NDIMZ];k++)
    for (j=0;j<lfmt[I_NDIMY];j++)
      for (i=0;i<lfmt[I_NDIMX];i++)
      {
         oldPhi = getValue(imgPhi,i,j,k);
         nvPhi = oldPhi + ptTau * getValue(imgF,i,j,k);
         setValue(imgPhi,i,j,k,nvPhi);
         if (nvPhi * oldPhi < 0)
           setValue(imgGeod,i,j,k,ptTau*(nb_iter - (nvPhi / (nvPhi - oldPhi))));
      }

  if (time_flag)
    {
      t1 = clock()-t0;
      printf("update_solution: temps mis: %f secondes\n",t1/(float)CLOCKS_PER_SEC);
    }
}

---