one_levelset.c

image processing including fourier,wavelet,segmentation etc.

		       if(old_dir=='h') {nb_angles++; break;}
		       if(old_dir=='b') {nb_angles--; break;}
		       break;}
	  }
	  len=lipptr->longueur;
	  while(--len>0){
	    newpc->x=(pX+=dX);
	    newpc->y=(pY+=dY);
	    newpc++;
	  }
	  pX+=dX;
	  pY+=dY;
	  old_dir=lipptr->direction;
	}                                            /* end-for lipptr     */
	if(old_dir!=first_dir)
	  switch(first_dir){
	    case 'b' : {if(old_dir=='g') {nb_angles++; break;}
			if(old_dir=='d') {nb_angles--; break;}
			break;}
	    case 'd' : {if(old_dir=='b') {nb_angles++; break;}
			if(old_dir=='h') {nb_angles--; break;}
			break;}
	    case 'h' : {if(old_dir=='d') {nb_angles++; break;}
			if(old_dir=='g') {nb_angles--; break;}
			break;}
	    case 'g' : {if(old_dir=='h') {nb_angles++; break;}
			if(old_dir=='b') {nb_angles--; break;}
			break;}
	  }
	assert((nb_angles==4)||(nb_angles==-4));
	/* as we have only closed Jordan curves the sum of angles in */
	/* a polygon is always 2\pi=4*(\pi/2), the sign gives the    */
	/* orientation of the curve (+=outer, - inner boundary)      */
	poly->channel[1]=(nb_angles==4)? (float)poly_nb : -(float)poly_nb ;
	poly++;/* this fixes the label&orientation  of the current polygon */
      }                                              /* end-for lbptr      */
      poly_nb++;
    }                                                /* end loop levelset  */
}



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

  for(rptr=image.tete;rptr!=NULL;rptr=rptr->Rsuiv) 
    if(*rptr->canal<=level) 
      for(lbptr=rptr->Lbords;lbptr!=NULL;lbptr=lbptr->suiv) 
        for (cnxeptr=lbptr->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);
}


void cDess_u_LS(boundary,u,level,inside,outside)  
float level;
char inside,outside;
Cimage boundary;          /* Balayage de l'image gauche/droite et haut/bas   */
Cimage 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]=(*regact->canal<=level)? inside:outside ;
                           /* 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]=(*regact->canal<=level)? inside:outside;
    }
  }
}


void fDess_u_LS(boundary,u,level)  
float level;
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-1;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]=*regact->canal;
                           /* Dessine le premier pixel de la i-ieme colonne   */
    for(j=1;j<nrows-1;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]=*regact->canal;
    }
  }
}


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
one_levelset(level,cb,pb,fu,cu,image_org)
float *level;
Cimage cb,cu;
Fpolygons pb;
Fimage fu;
Fimage image_org;
{
  Cimage boundary=NULL;
  unsigned long l;
  REGIONPTR rptr;
  LI_BORDSPTR lbptr;
  int number_lreg=0,number_points=0,number_cbords=0;

  if((!cb)&&(!pb)&&(!cu)&&(!fu))
    mwerror(USAGE,1,"Select at least one output option !\n");

  image.gx=image_org->ncol;
  image.gy=image_org->nrow;
  Estimation();
  Initialisation(image_org);

  segment(*level);

  for(rptr=image.tete;rptr!=NULL;rptr=rptr->Rsuiv)  
    if(*rptr->canal<=*level) number_lreg++;
  
  printf("\n Found %ld connected regions for the levelset.",number_lreg);

  if((cb)||(cu)||(fu)) {
    if((boundary=mw_change_cimage(NULL,image_org->ncol,image_org->nrow))==NULL)
      mwerror(FATAL,1,"Not enough memory.");
    for(l=0;l<image_org->nrow*image_org->ncol;l++) boundary->gray[l]=255;
    BlackBound(boundary,*level);
  }


  if(cu) {
    printf("\n Constructing b/w levelset picture.");
    cu=mw_change_cimage(cu,image_org->ncol,image_org->nrow);
    cDess_u_LS(boundary,cu,*level,0,255);
  }

  if(fu) {
    printf("\n Constructing gray (%.1f/%.1f) levelset picture.",
   MIN(*image.tete->Lbords->bord->rg->canal,*image.tete->Lbords->bord->rd->canal),
   MAX(*image.tete->Lbords->bord->rg->canal,*image.tete->Lbords->bord->rd->canal) );
    fu=mw_change_fimage(fu,image_org->ncol,image_org->nrow);
    fDess_u_LS(boundary,fu,*level);
  }

  if(cb) 
    *cb=*boundary;                               /* either spitt it out */
  else
    if((cu)||(fu)) mw_delete_cimage(boundary);   /* else to the trash   */
  
 
 if(pb) {
  /* Check that all regions have just connected, closed boundaries */
  /* else we will immediately exit, if everything is fine then     */
  /* count the total space needed in order to allocate memory (??) */
  /* Just consider the regions \leq than the current levellines    */
    for(rptr=image.tete;rptr!=NULL;rptr=rptr->Rsuiv)  
      if(*rptr->canal<=*level)
	for(lbptr=rptr->Lbords;lbptr!=NULL;lbptr=lbptr->suiv) {
	  if(lbptr->bord->cnxe->suiv!=NULL) {
	    fprintf(stderr,"\n Complex image structure... conversion to polygons not implemented\n");
	    Image_free(&image);
	    return;
	  }
	  number_points+=lbptr->bord->canal;
	  number_cbords++;
	}
  /* Now at least we are tranquil (hopefully!)                     */

    printf("\n Writting %ld contours, with %ld points, in MW2-polygons structure.",
	   number_cbords,number_points);
    pb=mw_change_fpolygons(pb);
    border2polygon(pb,number_cbords,*level);
  }
  Image_free(&image);
  printf("\n");
}