cml_decompose.c

image processing including fourier,wavelet,segmentation etc.

	    if(oldll!=NULL) oldll->previous=*lline;
	    newll->next=oldll;
	  }
      }

  /* open lines starting in the last row */
  for(c=0;c<NC-1;c++) 
    if((V[NL-1][c]==1)&&Inside(im,minvalue,maxvalue,NL-1,c)) 
      { /* ...levelset at left */
	oldll=*lline;
	/* Define a new lline but if newll corresponds to a former
	   computed lline with length < L
	*/
	if ((newll==NULL)||(newll==*lline) ) 
	  newll=produce_lline(minvalue,maxvalue,1);
	(*Nll)++;
	if ((p=newll->first_point)==NULL)
	  {
	    p=mw_new_point_curve();
	    if (p == NULL) mwerror(FATAL,1,"Not enough memory.\n");
	    newll->first_point=p;
	  }
	p->x=c+1;
	p->y=NL;
	V[NL-1][c]=0;
	sum=H[NL-2][c+1]+V[NL-2][c]+H[NL-2][c];
	follow_open_line(NL,NC,H,V,(int)NL-2,c,sum,newll);
	if (mw_length_cmorpho_line(newll) < L) (*Nllrm)++;
	else
	  {
	    *lline=newll;
	    if(oldll!=NULL) oldll->previous=*lline;
	    newll->next=oldll;
	  }
      }

  /* open lines starting in the first column */
  for(l=0;l<NL-1;l++)
    if((H[l][0]==1)&&Inside(im,minvalue,maxvalue,l,0)) 
      { /* ...levelset at left       */
	oldll=*lline;
	/* Define a new lline but if newll corresponds to a former
	   computed lline with length < L
	*/
	if ((newll==NULL)||(newll==*lline) ) 
	  newll=produce_lline(minvalue,maxvalue,1);
	(*Nll)++;
	if ((p=newll->first_point)==NULL)
	  {
	    p=mw_new_point_curve();
	    if (p == NULL) mwerror(FATAL,1,"Not enough memory.\n");
	    newll->first_point=p;
	  }
	p->x=0;
	p->y=l+1;
	H[l][0]=0;
	sum=V[l+1][0]+H[l][1]+V[l][0];
	follow_open_line(NL,NC,H,V,l,0,sum,newll);
	if (mw_length_cmorpho_line(newll) < L)  (*Nllrm)++;
	else
	  {
	    *lline=newll;
	    if(oldll!=NULL) oldll->previous=*lline;
	    newll->next=oldll;
	  }
      }
  
  /* open lines starting in the last column */
  for(l=0;l<NL-1;l++)
    if((H[l][NC-1]==1)&&!Inside(im,minvalue,maxvalue,l,NC-1)) 
      { /* ...levelset at left  */
	oldll=*lline;
	/* Define a new lline but if newll corresponds to a former
	   computed lline with length < L
	*/
	if ((newll==NULL)||(newll==*lline) ) 
	  newll=produce_lline(minvalue,maxvalue,1);
	(*Nll)++;
	if ((p=newll->first_point)==NULL)
	  {
	    p=mw_new_point_curve();
	    if (p == NULL) mwerror(FATAL,1,"Not enough memory.\n");
	    newll->first_point=p;
	  }
	p->x=NC;
	p->y=l+1;
	H[l][NC-1]=0;
	sum=V[l][NC-2]+H[l][NC-2]+V[l+1][NC-2];
	follow_open_line(NL,NC,H,V,l,(int)NC-2,sum,newll);
	if (mw_length_cmorpho_line(newll) < L) (*Nllrm)++;
	else
	  {
	    *lline=newll;
	    if(oldll!=NULL) oldll->previous=*lline;
	    newll->next=oldll;
	  }
    }
}

void
follow_closed_line(H,V,ll,cc,sum,lline)
unsigned char **H,**V;
int ll,cc,sum;
Cmorpho_line lline;
{
  Point_curve p0,p1;

  p0=lline->first_point;
  
  while(sum!=0) 
    {
      if ((p1=p0->next)==NULL)
	{
	  p1=mw_new_point_curve();
	  if (p1 == NULL) mwerror(FATAL,1,"Not enough memory.\n");
	  p1->previous=p0;
	  p0->next=p1;
	}
      p1->x=cc+1;
      p1->y=ll+1;
      p0=p1;
      if(sum==1) {
	if      (V[ll][cc]==1)     {V[ll][cc]=0;ll--;}
	else if (H[ll][cc]==1)     {H[ll][cc]=0;cc--;}
	else if (V[ll+1][cc]==1)   {V[ll+1][cc]=0;ll++;}
	else  /*(H[ll][cc+1]==1)*/ {H[ll][cc+1]=0;cc++;}
      }
      else { /* sum==3 */
	if     (V[ll][cc]==0)     {H[ll][cc+1]=0;cc++;}
	else if(H[ll][cc]==0)     {V[ll][cc]=0;ll--;}
	else if(V[ll+1][cc]==0)   {H[ll][cc]=0;cc--;}
	else /*(H[ll][cc+1]==0)*/ {V[ll+1][cc]=0;ll++;}
      }
      sum=V[ll][cc]+H[ll][cc]+V[ll+1][cc]+H[ll][cc+1];
  } /* end while */
}

void get_closed_lines(im,NL,NC,minvalue,maxvalue,H,V,lline,L,Nll,Nllrm)

Cfimage im;
Color minvalue,maxvalue;
unsigned int NL,NC;
unsigned char **H,**V;
Cmorpho_line *lline;
int L;
int *Nll,*Nllrm;

{
  Cmorpho_line oldll,newll=NULL;
  Point_curve p;
  int sum,l,c;

  for(l=0;l<NL-1;l++) for(c=0;c<NC-1;c++) 
    if(H[l][c+1]!=0) 
      {                    /* start a closed line */
	oldll=*lline;
	/* Define a new lline but if newll corresponds to a former
	   computed lline with length < L
	*/
	if ((newll==NULL)||(newll==*lline) ) 
	  newll=produce_lline(minvalue,maxvalue,0);
	(*Nll)++;
	if ((p=newll->first_point)==NULL)
	  {
	    p=mw_new_point_curve();
	    if (p == NULL) mwerror(FATAL,1,"Not enough memory.\n");
	    newll->first_point=p;
	  }
      p->x=c+1;
      p->y=l+1;
      if(!Inside(im,minvalue,maxvalue,l,c+1)) 
	{ /* such that the levelset is at left */
	  V[l+1][c]=0; 
	  sum=H[l+1][c]+V[l+2][c]+H[l+1][c+1];
	  follow_closed_line(H,V,l+1,c,sum,newll);
	}
      else 
	{
	  H[l][c+1]=0;
	  sum=V[l+1][c+1]+H[l][c+2]+V[l][c+1];
	  follow_closed_line(H,V,l,c+1,sum,newll);
	}
      if (mw_length_cmorpho_line(newll) < L) (*Nllrm)++;
      else
	{
	  *lline=newll;
	  if(oldll!=NULL) oldll->previous=*lline;
	  newll->next=oldll;
	}
      }
}

/******************************************************************************/
/*  Possible values for ml_opt :                                              */
/*  0  --> boundary{   levels[i]   <= im[x,y] <= MORPHO_INFTY } , level lines */
/*  1  --> boundary{-MORPHO_INFTY  <= im[x,y] <= levels[i]    } , inv. l. l.  */
/*  2  --> boundary{      levels[i]   == im[x,y]              } , iso lines   */
/*  3  --> boundary{   levels[i]   <= im[x,y] <= levels[i+1]  } , general ml  */
/*  4  --> boundary{  levels[2i]   <= im[x,y] <= levels[2i+1] } , id.         */
/*                                                       where i=0,1,...      */
/******************************************************************************/

void cml_extract(levels,opt,L,image_org,m_image)
Cfsignal levels;
int *opt;
int L;
Cfimage image_org;
Cmimage m_image;
{
  Cmorpho_line current_lline=NULL;
  unsigned long nb_points, l;
  unsigned long Nll, Nllrm;
  unsigned int NL=image_org->nrow, NC=image_org->ncol;
  unsigned char **V, **H, *cptr;
  int i,ml_opt=*opt;
  Color minvalue, maxvalue;

  /* memory for horizontal boundaries */
  H=(unsigned char **)malloc((NL-1)*sizeof(unsigned char*));
  H[0]=(unsigned char *)malloc((NL-1)*NC*sizeof(unsigned char));
  for(i=1;i<NL-1;i++) H[i]=H[i-1]+NC;

  /* memory for vertical boundaries */
  V=(unsigned char **)malloc(NL*sizeof(unsigned char*));
  V[0]=(unsigned char *)malloc(NL*(NC-1)*sizeof(unsigned char));
  for(i=1;i<NL;i++) V[i]=V[i-1]+(NC-1); 

  if((H[0]==NULL)||(V[0]==NULL)) {
    free((void*)H[0]);free((void*)V[0]);
    mwerror(FATAL,1,"Not enough memory.");
  }

  /* initialize V[][],H[][] to 0, i.e. no lines */
  cptr=V[0]+NL*(NC-1);      while(cptr-->*V) *cptr=0;
  cptr=H[0]+(NL-1)*NC;      while(cptr-->*H) *cptr=0;

  i=levels->size-1;
  Nll=Nllrm=0;
  do {
    if ((i%100)==0) mwdebug("--> Level i=%d\n",i);
    switch (ml_opt) 
      {
      case 0 : 
	minvalue=levels->values[i]; 
	maxvalue.red=maxvalue.green=maxvalue.blue=MORPHO_INFTY;  
	maxvalue.model=image_org->model;
	break;
      case 1 : 
	minvalue.red=minvalue.green=minvalue.blue=-MORPHO_INFTY; 
	minvalue.model=image_org->model;
	maxvalue=levels->values[i]; 
	break;
      case 2 : 
	minvalue=maxvalue=levels->values[i]; 
	break;
      default : mwerror(FATAL,1,"Bad option number");
      }

    produce_HV(image_org,NL,NC,minvalue,maxvalue,H,V);
    nb_points=count_X(H,V,NL,NC);
    /*
    printf("i=%d nb_points=%d L=%d\n",i,nb_points,L);
    */
    if(nb_points>=L) 
      {
	get_open_lines(image_org,NL,NC,minvalue,maxvalue,H,V,&current_lline,L,&Nll,&Nllrm);
	get_closed_lines(image_org,NL,NC,minvalue,maxvalue,H,V,&current_lline,L,&Nll,&Nllrm);
      }
    else   /* nb_points < L */
      {
	if (nb_points == 0)
	  {
	    switch (ml_opt) 
	      {
	      case 0 : case 1:
		mwerror(WARNING,1,
			"Value (%.5f,%.5f,%.5f) is smaller or bigger than values of image_in.\n",
			levels->values[i].red,levels->values[i].green,levels->values[i].blue);
		break;
		
	      case 2 : 
		mwerror(WARNING,1,"Value (%.5f,%.5f,%.5f) doesn't appear image_in.\n",
			levels->values[i].red,levels->values[i].green,levels->values[i].blue);
		break;
		
	      default : mwerror(FATAL,1,"Bad option number");
	      }
	    mwerror(WARNING,1,"No morpho_line generated.\n");
	  }
      }
    i--;
    if((ml_opt==3)&&(i==0))   i--;
    if((ml_opt==4)&&(--i==0)) i--;

  } while(i>=0);    

  mwdebug("Removed %d level lines of length < %d over %d (%3.1f %%).\n      Remain %d level lines.\n",Nllrm,L,Nll,(100.0*Nllrm)/Nll,Nll-Nllrm); 

  if(m_image==NULL) m_image=mw_change_cmimage(m_image);
  m_image->nrow=NL;
  m_image->ncol=NC;
  m_image->first_ml=current_lline;
  /* Other channels of mimage not to be completed by this module!! */

  free((void*)H[0]);free((void*)V[0]);
  free((void*)V);free((void*)H);
}

/* --------------------
   Main module function
   -------------------- */

Cmimage cml_decompose(cmimage_in,ml_opt,L,image_in)
Cmimage cmimage_in;
int* ml_opt;
int *L;
Cfimage image_in;
{
  Cmimage cmimage=NULL;
  Cfsignal levels,tmp_levels;
  long l;
  
  cmimage=mw_change_cmimage(cmimage);
  if (cmimage == NULL) mwerror(FATAL,1,"Not enough memory.\n");

  if(cmimage_in) 
    {
      if(cmimage_in->first_ml!=NULL)
	mwerror(WARNING,1,"Level lines of cmimage_in not copied!\n");
      if((cmimage_in->nrow!=image_in->nrow)||
	 (cmimage_in->ncol!=image_in->ncol))
	mwerror(WARNING,1,"image_in and cmimage_in not of same dimension!\n Dimensions of image_in are kept.\n");
      cmimage->first_fml=cmimage_in->first_fml;
      cmimage->first_ms=cmimage_in->first_ms;
    }
  else 
    {
      cmimage->first_fml=NULL;
      cmimage->first_ms=NULL;
    }
  cmimage->nrow=image_in->nrow;
  cmimage->ncol=image_in->ncol;

  /* Color of input image should be in MODEL_HSI */
  if (image_in->model != MODEL_HSI)
    mwerror(WARNING,0,"Color model of input cfimage is not HSI !\n");

  /* construct levels according to the ml_opt */
  if(*ml_opt==1) ascending_order=-1; else ascending_order=1;
  levels = cfvalues(image_in);
  cmimage->minvalue=levels->values[0];
  cmimage->maxvalue=levels->values[levels->size-1];
  if (*ml_opt!=2)
    {
      levels->size--;
      for(l=0;l<levels->size;l++) levels->values[l]=levels->values[l+1];
    }

  if(2*levels->size>cmimage->nrow*cmimage->ncol)
    printf("\n Warning : %d different pixel values (for %d pixels) !\n",levels->size,cmimage->nrow*cmimage->ncol);

  cml_extract(levels,ml_opt,*L,image_in,cmimage);

  free(levels->values);
  free(levels);

  if (mwdbg == 1)
    {
      mwdebug("Checking mimage in ml_decompose (%d level lines)...\n",
              mw_num_cmorpho_line(cmimage->first_ml));
      llcheck(cmimage);
      mwdebug("End of checking mimage\n");
    }
  return(cmimage);
}