one_levelset.c

image processing including fourier,wavelet,segmentation etc.

  for (li=bord1->rg->Lbords;li->bord!=bord1;li=li->suiv) ;
  if (       ((li->suiv!=NULL) && (li->suiv->bord!=bord2)) 
        ||   ((li->suiv==NULL) && (bord2!=bord1->rg->Lbords->bord))   )
    {
         bord =bord1;       bord1=bord2;
         bord2=bord ;
    }


  /* Determination des comp. cnxes. de sommet 'som' et de leur orientation */
  /* avec:                                                                 */
  /*                                                              */
  /*         bcnxe1            som        bcnxe2         */
  /*     sb -------> sa (+1)   $     sb -------> sa (+1) */
  /*     sb -------> sa (+1)   $     sa <------- sb (-1) */
  /*     sa <------- sb (-1)   $     sb -------> sa (+1) */
  /*     sa <------- sb (-1)   $     sa <------- sb (-1) */
  /*                                                     */
  /*******************************************************/

  for(bcnxe=bord1->cnxe;bcnxe!=NULL;bcnxe=bcnxe->suiv) {
     if (bcnxe->sa==som) {bcnxe1=bcnxe;dir_bcnxe1=1;break;}
     if (bcnxe->sb==som) {bcnxe1=bcnxe;dir_bcnxe1=-1;break;}
  }

  for(bcnxe=bord2->cnxe;bcnxe!=NULL;bcnxe=bcnxe->suiv) {
     if (bcnxe->sa==som) {bcnxe2=bcnxe;dir_bcnxe2=-1;break;}
     if (bcnxe->sb==som) {bcnxe2=bcnxe;dir_bcnxe2=1;break;}
  }

 /* Mettre dans bord1 les deux bords: b1 <- b1+b2 */

  bord1->canal += bord2->canal ;

  /* Est-ce qu'il faut changer l'orientation des bcnxes de bord2 */
  if (dir_bcnxe1!=dir_bcnxe2) 
      for(bcnxe=bord2->cnxe;bcnxe!=NULL;bcnxe=bcnxe->suiv) {
          sa_old=bcnxe->sa; bcnxe->sa=bcnxe->sb; bcnxe->sb=sa_old;
          lab_old=bcnxe->lab;bcnxe->lab=bcnxe->lba;bcnxe->lba=lab_old;
      }

  /* Repointer les sommets sur bord */
  for(bcnxe=bord2->cnxe;bcnxe!=NULL;bcnxe=bcnxe->suiv) {
    Repointer(bcnxe->sa,bord1,bord2);
    Repointer(bcnxe->sb,bord1,bord2);
  }

  /* Mettre les bcnxes ensembles dans bord1 */
  for(bcnxe=bord1->cnxe;bcnxe->suiv!=NULL;bcnxe=bcnxe->suiv) ;
  bcnxe->suiv=bord2->cnxe;

  /* Reunir les deux bcnxes dont som est un sommet dans bcnxe1 */
  if (dir_bcnxe1 == 1) {
         bcnxe1->sa=bcnxe2->sa;
         bcnxe1->lba=LiPixelsUnion(bcnxe1->lba,bcnxe2->lba);
         bcnxe1->lab=LiPixelsUnion(bcnxe2->lab,bcnxe1->lab);
  }
  else {
         bcnxe1->sb=bcnxe2->sb;
         bcnxe1->lba=LiPixelsUnion(bcnxe2->lba,bcnxe1->lba);
         bcnxe1->lab=LiPixelsUnion(bcnxe1->lab,bcnxe2->lab);
  }
  for(bcnxe=bord1->cnxe;bcnxe->suiv!=bcnxe2;bcnxe=bcnxe->suiv) ;
  bcnxe->suiv=bcnxe2->suiv;

  /* Elimination de bord2 des listes de bords des regions */
  ElimBordeReg(bord1->rg,bord2);
  ElimBordeReg(bord1->rd,bord2);
}

void UnionBordCnxe(reg)  /* Reunit les bords de reg,qui font frontiere a une */
REGIONPTR reg;           /*  meme region voisine,                            */
{                        /*  on obtient le comp. connexes d'un bord          */
  short i,dir_bord2;
  static void Repointer();
  SOMMETPTR som;
  LI_PIXELSPTR lpix;
  LI_BORDSPTR lbor1,lbor2,lbor3,lbor;
  BORDPTR bord1,bord2;
  BORDCONNEXEPTR bcnxe;
  REGIONPTR reg1,reg2;
  
  lbor1=reg->Lbords;
  while((lbor1!=NULL)&&(lbor1->suiv!=NULL)) {
    bord1=lbor1->bord;
    reg1 = (bord1->rg==reg) ?  bord1->rd : bord1->rg ;
    lbor2=lbor1->suiv;
    do {
      bord2=lbor2->bord;
      reg2 = (bord2->rg==reg) ?  bord2->rd : bord2->rg ;
     /*Est-ce qu'il y a deux bords de reg,frontiere avec une meme autre region*/
      if (reg1==reg2) {
         image.nbbords--;                     /* De deux bords on en fait un */
         dir_bord2 = (bord1->rg==bord2->rg) ?  1 : -1 ;
         /* Les sommets des comp. cnxes. de bord2 doivent pointer vers bord1 */
         for(bcnxe=bord2->cnxe;bcnxe!=NULL;bcnxe=bcnxe->suiv) {
            Repointer(bcnxe->sa,bord1,bord2);
            Repointer(bcnxe->sb,bord1,bord2);
            if (dir_bord2==-1) {
               som=bcnxe->sa;bcnxe->sa=bcnxe->sb;bcnxe->sb=som;
               lpix=bcnxe->lab;bcnxe->lab=bcnxe->lba;bcnxe->lba=lpix;
            }  
         }
         /* Mettre les parties connexes dans bord1 */
         bord1->canal+=bord2->canal;
	 for(bcnxe=bord1->cnxe;bcnxe->suiv!=NULL;bcnxe=bcnxe->suiv) ;
         bcnxe->suiv=bord2->cnxe;

         /* Eliminer bord2 de la liste de bords de reg */
         for(lbor3=lbor1;lbor3->suiv->bord!=bord2;lbor3=lbor3->suiv) ;
         lbor=lbor3->suiv;
         lbor3->suiv=lbor->suiv;

         /* Eliminer bord2 de la liste de bords de reg1=reg2 */
         lbor3=reg1->Lbords;
         if (lbor3->bord==bord2) {
             lbor=lbor3;
             reg1->Lbords=lbor->suiv;
         }
         else {
             for( ;lbor3->suiv->bord!=bord2;lbor3=lbor3->suiv) ;
             lbor=lbor3->suiv;
             lbor3->suiv=lbor->suiv;
         }
         break;
      }
      else lbor2=lbor2->suiv;
    } while(lbor2!=NULL);
    lbor1=lbor1->suiv;
  }    
}


void Repointer(som,newbord,oldbord) /* som va pointer sur newbord et plus sur*/
SOMMETPTR som;                      /*  oldbord,si'newbord'est deja dans la  */
BORDPTR newbord,oldbord;            /*  liste som->Lbords on elimine'oldbord'*/
{
  LI_BORDSPTR li;
  unsigned char existe_new;

  existe_new=FAUX;
  for(li=som->Lbords;li!=NULL;li=li->suiv) 
            existe_new |=(li->bord==newbord);
  if (existe_new==VRAI) {  /* On elimine 'oldbord' de la liste 'som->Lbords' */
       li=som->Lbords;
       if(li->bord==oldbord) 
         som->Lbords=li->suiv;
       else {
         for(li=som->Lbords;(li->suiv!=NULL)&&(li->suiv->bord!=oldbord);li=li->suiv) ;
         if(li->suiv==NULL) return;  /*'newbord existe deja et 'oldbord' plus */
         li->suiv=li->suiv->suiv;
       }
  }
  else  {               /* On fait pointer 'som' sur 'newbord' */
       for(li=som->Lbords;li->bord!=oldbord;li=li->suiv) ;
       li->bord=newbord;
  } 
}


void ElimBordeReg(reg,bord)                /* Elimination de bord dans la */
REGIONPTR reg;                             /*  liste de bords de region   */
BORDPTR bord;
{
  LI_BORDSPTR li;
 
  li=reg->Lbords;
  if (li->bord==bord)
     reg->Lbords=li->suiv;
  else {
     for ( ;li->suiv->bord!=bord;li=li->suiv) ;
         li->suiv=li->suiv->suiv;
     }
}



LI_PIXELSPTR 
LiPixelsUnion(liste1,liste2)                /* Met la liste de pixels liste2 */
LI_PIXELSPTR liste1,liste2;                 /*  derriere liste1              */
{
  LI_PIXELSPTR liste,lptr;

  for(liste=liste1;liste->suiv!=NULL;liste=liste->suiv);
  if (liste->direction==liste2->direction) {
    liste->longueur += liste2->longueur;
    lptr=liste2;
    liste->suiv=lptr->suiv;
  }
  else liste->suiv=liste2;

  return(liste1);
}


/* 'segmentation' PROCEDURES */
void
segment(level)
float level;
{
  static void RegMerge();
  unsigned char agglutination=VRAI;
  REGIONPTR regact,regvois;
  BORDPTR bordcom;
  LI_BORDSPTR libords;

  /* Boucle balayage de l'image jusqu'a obtention des composantes connexes */
  while (agglutination==VRAI){
    agglutination=FAUX;
    for (regact=image.tete;regact!=NULL;regact=regact->Rsuiv) {
      for (libords=regact->Lbords;libords!=NULL;libords=libords->suiv) {
        bordcom=libords->bord;
        regvois=(bordcom->rd==regact)   ?   bordcom->rg : bordcom->rd ;
        if(*regact->canal==*regvois->canal) {
          agglutination=VRAI;
          RegMerge(regact,regvois,bordcom);
          break;                               /* we found a candidate region */
        }
      }
    }
  }
  printf(" Decomposition in connected components finished.\n");
  agglutination=VRAI;

  printf(" Extraction of contours of regions with gray <= %.4f.\n",level);
  while (agglutination==VRAI){
    agglutination=FAUX;
    for (regact=image.tete;regact!=NULL;regact=regact->Rsuiv) {
      for (libords=regact->Lbords;libords!=NULL;libords=libords->suiv) {
        bordcom=libords->bord;
        regvois=(bordcom->rd==regact)   ?   bordcom->rg : bordcom->rd ;
        if(  ((*regact->canal<=level)&&(*regvois->canal<=level))
           ||((*regact->canal> level)&&(*regvois->canal> level)) ){
          agglutination=VRAI;
          RegMerge(regact,regvois,bordcom);
          break;                             /* we found a candidate region */
        }
      }
    }
  }
}

/* This function puts the coordinates of the points from the segmentation    */
/* in 'point_curves'. We give the (integer) coordinates of the tips...       */
void
border2polygon(setpoly,nb_poly,level) 
Fpolygons setpoly;
int nb_poly;
float level;
{
  float pX,pY,dX,dY;
  int i,nb_points,nb_angles;
  unsigned short len,poly_nb=1;
  float *f;
  REGIONPTR rptr;
  BORDPTR bord;
  LI_BORDSPTR lbptr;
  LI_PIXELSPTR lipptr;
  Fpolygon poly;
  Point_fcurve newpc;
  char old_dir,first_dir;

  if((setpoly->first=(Fpolygon)malloc(nb_poly*sizeof(struct fpolygon)))==NULL)
    mwerror(FATAL,1,"Not enough memory.");
  if((f=(float*)malloc(2*nb_poly*sizeof(float)))==NULL)
    mwerror(FATAL,1,"Not enough memory.");

  /* Initialize the 'nb_poly' polygons */
  for(i=0,poly=setpoly->first;i<nb_poly;i++,poly++) {
    poly->nb_channels=2;
    poly->channel=f++;f++;
    poly->first=NULL;
    poly->previous=(i==0)? NULL:(poly-1) ;
    poly->next=(i==nb_poly-1)? NULL:(poly+1) ;
  }

  poly=setpoly->first;
  for(rptr=image.tete;rptr!=NULL;rptr=rptr->Rsuiv) 
    if(*rptr->canal<=level) {
      for(lbptr=rptr->Lbords;lbptr!=NULL;lbptr=lbptr->suiv) {
	bord=lbptr->bord; 
	poly->channel[0]=level;    /* here we give the levelset to the polygon */
	nb_points=(int)bord->canal;
	if((poly->first=(Point_fcurve)
	           malloc(nb_points*sizeof(struct point_fcurve)))==NULL)
	  mwerror(FATAL,1,"Not enough memory.");
	/* Initialize the 'nb_points' points */
	for(i=0,newpc=poly->first;i<nb_points;i++,newpc++) {
	  newpc->previous=(i==0)? NULL:(newpc-1) ;
	  newpc->next=(i==nb_points-1)? NULL:(newpc+1) ;
	}
	if(bord->rg==rptr) {
	  pX=(float)bord->cnxe->sb->x-0.5;
	  pY=(float)bord->cnxe->sb->y-0.5;
	  lipptr=bord->cnxe->lba;
	}
	else {
	  pX=(float)bord->cnxe->sa->x-0.5;
	  pY=(float)bord->cnxe->sa->y-0.5;
	  lipptr=bord->cnxe->lab;
	}
	first_dir=lipptr->direction;
	nb_angles=0;
	old_dir='0';
	newpc=poly->first;
	for ( ;lipptr!=NULL;lipptr=lipptr->suiv){
	  newpc->x=pX;
	  newpc->y=pY;
	  newpc++;
	  switch(lipptr->direction){
	   case 'b' : {dX=  0.0;dY=  1.0 ;
		       if(old_dir=='g') {nb_angles++; break;}
		       if(old_dir=='d') {nb_angles--; break;}
		       break;}
	   case 'd' : {dX=  1.0;dY=  0.0 ;
		       if(old_dir=='b') {nb_angles++; break;}
		       if(old_dir=='h') {nb_angles--; break;}
		       break;}
	   case 'h' : {dX=  0.0;dY= -1.0 ;
		       if(old_dir=='d') {nb_angles++; break;}
		       if(old_dir=='g') {nb_angles--; break;}
		       break;}
	   case 'g' : {dX= -1.0;dY=  0.0 ;