segct.c

image processing including fourier,wavelet,segmentation etc.

  rg=bord->rg;
  rd=bord->rd;
  coef_aire =rg->canal[0]*rd->canal[0]/(rg->canal[0]+rd->canal[0]);
  u1_2=(rg->canal[1]/rg->canal[0])-(rd->canal[1]/rd->canal[0]);
  som_p += u1_2 * u1_2 ;
  t_elastique =  coef_aire * som_p ;
  return( t_elastique / bord->canal[0] );
}





/* 'segmentation' PROCEDURES */


void
segment()
{
  float lam;
  BORDPTR bord;
  static void RegMerge();

  lam=h_lambda(h_root);   /* get the lowest lambda value */
  image.lambda=lam;
  do{
    bord=h_bord(h_root);/* get the lightest border */
    image.energie.l-=(long)bord->canal[0];
    image.energie.e+=(float)bord->canal[0] * h_lambda(h_root);
    RegMerge(bord->rg,bord->rd,bord);
  }while(image.nbbords>0 && lam>=h_lambda(h_root));
  return;
}



void EcritComm(comm)     /* Ecrit le commentaire pour le fichier resultat,   */
char *comm;              /*  ceci suivant les particularites de l'algorithme.*/
{                        /* Utilisee par 'SaveIm()' et autres 'Save'         */

   sprintf(comm,"2-normal segmentation (Mumford&Shah model), grey-level. Parameter:%lu, grid of %d*%d pixels. Nb. of regions: %d,frontiers: %d,tips: %d, El.en. %.4g,length %ld(G.Koepfler).",image.lambda,image.PAS,image.PAS,image.nbregions,image.nbbords,image.nbsommets,image.energie.e,image.energie.l);
}


make_curves(curves)
Curves curves;
{
  REGIONPTR regptr;
  LI_BORDSPTR libptr;
  BORDCONNEXEPTR cnxeptr;
  LI_PIXELSPTR lipptr;
  Curve oldcv,newcv;
  Point_curve oldpc,newpc;
  short a0,b0,a1,b1,dx,dy;
  
  newcv = oldcv = NULL;

  for (regptr=image.tete; regptr!=NULL; regptr=regptr->Rsuiv)
    for (libptr=regptr->Lbords;libptr!=NULL;libptr=libptr->suiv)
      for (cnxeptr=libptr->bord->cnxe;cnxeptr!=NULL;cnxeptr=cnxeptr->suiv)
	{
	  oldcv = newcv;	  /* A new curve */
	  newcv = mw_new_curve();
	  if (newcv == NULL)
	    {
	      mw_delete_curves(curves);
	      mwerror(FATAL,1,"Not enough memory\n");
	    }
	  if (curves->first == NULL) curves->first = newcv;
	  if (oldcv != NULL) oldcv->next = newcv;
	  newcv->previous = oldcv;
	  newcv->next = NULL;

	  oldpc = newpc = NULL;
	  newpc = mw_new_point_curve();
	  if (newpc == NULL)
	    {
	      mw_delete_curves(curves);
	      mwerror(FATAL,1,"Not enough memory\n");
	    }
	  if (newcv->first == NULL) newcv->first = newpc;
	  if (oldpc != NULL) oldpc->next = newpc;
	  newpc->previous = oldpc;
	  newpc->next = NULL;	  
	  newpc->x = a0 = cnxeptr->sb->x;
	  newpc->y = b0 = cnxeptr->sb->y;

	  for (lipptr=cnxeptr->lba;lipptr!=NULL;lipptr=lipptr->suiv)
	    {
	      dx=dy=0;
	      switch (lipptr->direction)
		{
		case 'b' : dy=  lipptr->longueur ;break;
		case 'd' : dx=  lipptr->longueur ;break;
		case 'h' : dy= -lipptr->longueur ;break;
		case 'g' : dx= -lipptr->longueur ;break;
		}
	      a1=a0+dx;
	      b1=b0+dy;
	      
	      oldpc = newpc;
	      newpc = mw_new_point_curve();
	      if (newpc == NULL)
		{
		  mw_delete_curves(curves);
		  mwerror(FATAL,1,"Not enough memory\n");
		}
	      if (newcv->first == NULL) newcv->first = newpc;
	      if (oldpc != NULL) oldpc->next = newpc;
	      newpc->previous = oldpc;
	      newpc->next = NULL;	  
	      newpc->x = a0 = a1;
	      newpc->y = b0 = b1;
	    }
	}	      
}


void Dess_u(boundary,u)  
Cimage boundary;          /* Balayage de l'image gauche/droite et haut/bas   */
Fimage u;                 /* Determine s'il ya un bord ou non                */
{                         /* Si oui on change de region, dans les deux cas   */
  REGIONPTR regact,regtop;/*on dessine en utilisant la fonction correspondant*/
  static REGIONPTR RegAdjD(),RegAdjB();           /* a la region actuelle.   */
  short i,j,nrows=boundary->nrow,ncols=boundary->ncol;

  regtop=image.regpixel0_0;
  for(i=0;i<ncols;i++) {
    regact=regtop;
    if((i!=ncols-1)&&(boundary->gray[i+1]==0)) 
      regtop=RegAdjD(regact,i+1,0); /* Fin image ou bord vertical 1ere ligne */
    u->gray[i]=(float)regact->canal[1]/regact->canal[0];
                           /* Dessine le premier pixel de la i-ieme colonne   */
    for(j=1;j<nrows;j++){/* Dessine les pixels restants de la i-ieme colonne*/
      if(boundary->gray[ncols*j+i]==0)  regact=RegAdjB(regact,i,j) ;
      u->gray[ncols*j+i]=(float)regact->canal[1]/regact->canal[0];
    }
  }
}


REGIONPTR RegAdjD(reg,i,j)
REGIONPTR reg;
short i,j;
{
  LI_BORDSPTR lbords;
  BORDPTR bordptr;
  BORDCONNEXEPTR cnxeptr;
  LI_PIXELSPTR lipptr;
  short a0,b0,a1,b1,dx,dy;

  for(lbords=reg->Lbords;lbords!=NULL;lbords=lbords->suiv)  {
    bordptr=lbords->bord;
    cnxeptr=bordptr->cnxe;
    while(cnxeptr!=NULL)  {
      a0=cnxeptr->sb->x;
      b0=cnxeptr->sb->y;
      for (lipptr=cnxeptr->lba; lipptr!=NULL; lipptr=lipptr->suiv) {
        dx=dy=0;
        switch (lipptr->direction){
        case 'b' : dy=   lipptr->longueur; break;
        case 'd' : dx=   lipptr->longueur; break;
        case 'h' : dy=  -lipptr->longueur; break;
        case 'g' : dx=  -lipptr->longueur; break;
        }
        a1=a0+dx;
        b1=b0+dy; /* conditions pour que (i,j) soit un pixel de bord vertical */
        if ((a0==a1)&&(i==a0)&&( ((j>=b1)&&(j<b0))||((j>=b0)&&(j<b1)) ))
          return((reg==bordptr->rg) ? bordptr->rd : bordptr->rg ) ;
        a0=a1;b0=b1;
      }
      cnxeptr=cnxeptr->suiv;
    }
  }
  return(reg);
}


REGIONPTR RegAdjB(reg,i,j)
REGIONPTR reg;
short i,j;
{
  LI_BORDSPTR lbords;
  BORDPTR bordptr;
  BORDCONNEXEPTR cnxeptr;
  LI_PIXELSPTR lipptr;
  short a0,b0,a1,b1,dx,dy;
  
  for(lbords=reg->Lbords;lbords!=NULL;lbords=lbords->suiv)  {
    bordptr=lbords->bord;
    cnxeptr=bordptr->cnxe;
    while(cnxeptr!=NULL)  {
      a0=cnxeptr->sb->x;
      b0=cnxeptr->sb->y;
      for (lipptr=cnxeptr->lba; lipptr!=NULL; lipptr=lipptr->suiv) {
        dx=dy=0;
        switch (lipptr->direction){
        case 'b' : dy=   lipptr->longueur; break;
        case 'd' : dx=   lipptr->longueur; break;
        case 'h' : dy=  -lipptr->longueur; break;
        case 'g' : dx=  -lipptr->longueur; break;
        }
        a1=a0+dx;
        b1=b0+dy;/* conditions pour que (i,j) soit un pixel de bord horizontal*/
        if ((b0==b1)&&(j==b0)&&( ((i>=a1)&&(i<a0))||((i>=a0)&&(i<a1)) ))
          return((reg==bordptr->rg) ? bordptr->rd : bordptr->rg ) ;
        a0=a1;b0=b1;
      }
      cnxeptr=cnxeptr->suiv;
    }
  }
  return(reg);
}

 
void BlackBound(whitesheet)
Cimage whitesheet;
{
  static short TraitHVmono();
  short a0,b0,a1,b1,dx,dy;
  REGIONPTR regptr;
  LI_BORDSPTR libptr;
  BORDCONNEXEPTR cnxeptr;
  LI_PIXELSPTR lipptr;

  printf("\n Writting the boundaries.\n");
  for (regptr=image.tete; regptr!=NULL; regptr=regptr->Rsuiv){
    for (libptr=regptr->Lbords;libptr!=NULL;libptr=libptr->suiv){
      if (libptr->bord->rg!=regptr){
        for (cnxeptr=libptr->bord->cnxe;cnxeptr!=NULL;cnxeptr=cnxeptr->suiv){
          a0=cnxeptr->sb->x;
          b0=cnxeptr->sb->y;
          for (lipptr=cnxeptr->lba;lipptr!=NULL;lipptr=lipptr->suiv){
            dx=dy=0;
            switch (lipptr->direction){
             case 'b' : dy=  lipptr->longueur ;break;
             case 'd' : dx=  lipptr->longueur ;break;
             case 'h' : dy= -lipptr->longueur ;break;
             case 'g' : dx= -lipptr->longueur ;break;
            }
            a1=a0+dx;
            b1=b0+dy;
            TraitHVmono(whitesheet,a0,b0,a1,b1,0);
            a0=a1;
            b0=b1;
          }
        }
      }
    }
  }
}


short TraitHVmono(whitesheet,a0,b0,a1,b1,c)
Cimage whitesheet;
short a0, b0,a1,b1;
unsigned char c;
{
  short bdx=whitesheet->ncol,
        bdy=whitesheet->nrow;
  short z=0,dz,sz;

  if ((a0!=a1)&&(b0!=b1)) return(-2);
  if (b0==b1) {
    if (a0<0) a0=0; else if (a0>=bdx) a0=bdx-1;
    if (a1<0) {a1=0; if (a0==0) return(-1);}
    else if (a1>=bdx) {a1=bdx-1; if (a0==bdx-1) return(-1);}
    if (a0<a1) {sz=1; dz=a1-a0; }
    else {sz=-1; dz=a0-a1; }
    while (abs(z)<=dz) {
      whitesheet->gray[((long) bdx)*b0+z+a0]=c;
      z+=sz;
    }
  }
  else {
    if (b0<0) b0=0; else if (b0>=bdy) b0=bdy-1;
    if (b1<0) { b1=0; if (b0==0) return(-1); }
    else if (b1>=bdy) {b1=bdy-1; if (b0==bdy-1) return(-1); }
    if (b0<b1) {sz=1; dz=b1-b0; } else {sz=-1; dz=b0-b1; }
    while (abs(z)<=dz) {
      whitesheet->gray[((long) bdx)*(z+b0)+a0]=c;
      z+=sz;
    }
  }
  return(0);
}


/* Main function */
 
Cimage segct(sgrid,nb_of_regions,lambda,curves,u,f_nb_of_regions,f_lambda,image_org)
int *sgrid,*nb_of_regions,*f_nb_of_regions;
float *lambda,*f_lambda;
Fimage image_org,u;
Curves curves;
{
  Cimage boundary;
  unsigned long l;

  if((nb_of_regions==NULL)==(lambda==NULL))
    mwerror(FATAL,1,"Chose exactly one option out of '-N' and '-L'.\n");

  if(sgrid) image.PAS = *sgrid;
  else      image.PAS = 1;
  image.NC_REG=2;
  /*
  if((image_org->ncol % image.PAS != 0) || (image_org->nrow % image.PAS != 0)) 
  */
    if ((image_org->ncol != (image.PAS*(image_org->ncol/image.PAS))) || 
        (image_org->nrow != (image.PAS*(image_org->nrow/image.PAS))) )
      mwerror(FATAL,1,"Use an image with dimensions multiples of %d !\n",image.PAS);

  image.gx = image_org->ncol / image.PAS;
  image.gy = image_org->nrow / image.PAS;

  Estimation();
  Initialisation(image_org);
  printf("\nInitial state of the segmentation:");
  printf("\nElastic energie = %.4g , boundary length= %ld.",image.energie.e,image.energie.l);
  printf("\nNumber of regions: %d. \n",image.nbregions);
  if(lambda==NULL) 
    while(*nb_of_regions<image.nbregions) segment();
  else 
    while(*lambda>=h_lambda(h_root)) segment(); /* h_lambda(h_root) contains   */
                                      /* the new value for which we will merge */
  *f_lambda=image.lambda;
  *f_nb_of_regions=image.nbregions;
  printf("\nElastic energie = %.4g , boundary length= %ld.",image.energie.e,image.energie.l);
  printf("\nNumber of regions: %ld, value of lambda %f \n",image.nbregions,image.lambda);

  if((boundary=mw_change_cimage(NULL,image_org->nrow,image_org->ncol))==NULL)
    {Image_free(&image);mwerror(FATAL,1,"Not enough memory.");}
  for(l=0;l<image_org->nrow*image_org->ncol;l++) boundary->gray[l]=255;
  BlackBound(boundary);
  if(u) {
    if((u=mw_change_fimage(u,image_org->nrow,image_org->ncol))==NULL)
      {Image_free(&image);mwerror(FATAL,1,"Not enough memory !\n");}
    Dess_u(boundary,u);
  }
  if(curves) make_curves(curves);
  Image_free(&image);
  return(boundary);
}