ll_boundaries2.c

image processing including fourier,wavelet,segmentation etc.

/*----------------------------- MegaWave Module -----------------------------*/
/* mwcommand
   name = {ll_boundaries2};
   version = {"1.0"};
   author = {"Lionel Moisan, Frederic Cao"};
   function = {"Extract (local or not) meaningful boundaries from a Fimage"};
   usage = {
'e':[eps=0.]->eps   "-log10(max number of false alarms), default 0.",
's':[step=1.]->step "quantization step, default 1.",
'p':[prec=2]->prec  "sampling precision for flst_bilinear, default 2",
'G':[G=0.5]->std    "standard dev. for preliminary convolution, default 0.5",
'a'->all            "keep all meaningful level lines (not only maximal ones)",
'H':[hstep=0.01]->hstep "step for contrast histogram, default 0.01",
't':tree->tree          "use a precomputed bilinear tree, default NULL", 
'v':[visit=100]->visit  "maximal number of visits for a boundary, default 100",
'L'->loc                  "force local research",
'k':keep_tree<-keep_tree  "keep meaningful tree",
'o':image_out<-image_out  "image reconstructed by meaningful tree",
in->in                    "input (Fimage)",
out<-ll_boundaries2       "output boundaries (Flists)"
}
;
*/

#include <stdio.h>
#include <float.h>
#include <math.h>
#include "mw.h"

extern void flst();
extern void flst_bilinear();
extern void flst_reconstruct();
extern void flstb_quantize();
extern Flist flstb_boundary();
extern void flst_pixels();
extern int  fsaddles();
extern Fsignal sintegral();
extern Fimage fsepconvol();

#define ABS(x) ((x)>0?(x):-(x))
#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MIN_HIST_AREA 10000

typedef struct mydata{
  int visit;
  int ndetect; 
  char type;
  float nfa; /* number of false alarms */
  float min_contrast; 
  float p_contrast;
  float length;
  Fsignal histo; /* histogram of gradient inside the shape */
  char meaningful; 
  char max; /* maximal meaningful boundary */
  char notmax; /* meaningful boundary known to be NOT maximal */
} *Mydata;

/* A type containing fixed values. This prevents from global variables 
   and huge function declaration.   */


struct myglobal{
  Shapes Tree;
  float **tabsaddles,threshold;
  Fimage image;
  int *prec;
  char *local,*all;
};



/* store meaningful boundaries */
void store_boundaries(cs,tree,image,prec,tabsaddles,visit,all,
		      NormofDu,eps2)
Flists cs;
Shapes tree;
int *prec,visit;
float **tabsaddles,eps2;
Fimage image,NormofDu;
char *all;
{
  Shape s;
  Flist l,lhisto;
  int i,j;
  Mydata sdata; 
  float proba,probastep,lognbtests,proba2;
  

  lognbtests = (float)log10((double)tree->nb_shapes);
  proba2 = (float) pow(10.,(double)eps2);
  l = mw_new_flist();
  if(!l) mwerror(FATAL,1,"Not enough memory.\n");
  for(i=1;i<tree->nb_shapes;i++){
    s = tree->the_shapes+i;
    sdata = (Mydata) s->data;
    if(sdata->max || (all && sdata->meaningful)){   
      s->removed = 0;
      mwdebug("#%d shape #%d c=%.2f pc=%.4f l=%.2f PFA=%.3f\n",
	      cs->size,i,sdata->min_contrast,sdata->p_contrast,
	      sdata->length,sdata->nfa);  
      flstb_boundary(prec,image,tree,s,NULL,l,tabsaddles);
      cs->list[cs->size] = mw_copy_flist(l,NULL);
      
      /* store nfa as a double in the data field */
      cs->list[cs->size]->data_size = sizeof(double);
      cs->list[cs->size]->data = malloc(sizeof(double));
      *((double*)(cs->list[cs->size]->data)) = 
	(double)(sdata->nfa+lognbtests);
      cs->size++;
    }
    else s->removed = 1;
  }
  mw_delete_flist(l);
  
  /* store the initial number of curve in the data field of the Flists*/
  cs->data_size = sizeof(double);
  cs->data = malloc(sizeof(double));
  *((double*)cs->data) = (double) tree->nb_shapes;
}


/*===== Compute the minimum contrast and the length of the curve l =====*/

float min_contrast(l,length,NormofDu)
Flist l;
float *length;
Fimage NormofDu;
{
  double per;
  float mu,minmu,x,y,ox,oy;
  int i,ix,iy,minx,miny;

  per = 0.;
  minmu = FLT_MAX;

  for (i=0;i<l->size;i++) {

    x = l->values[i*2];
    y = l->values[i*2+1];

    if (i>0) per += sqrt((double)(x-ox)*(x-ox)+(y-oy)*(y-oy));
    ox = x; oy = y;
    
    ix = (int)rint((double)x)-1;
    iy = (int)rint((double)y)-1;
    if (ix>=0 && iy>=0 && ix<NormofDu->ncol && iy<NormofDu->nrow) {
      mu = NormofDu->gray[NormofDu->ncol*iy+ix];
      if (mu<minmu) {
	minmu= mu;
	minx = ix;
	miny = iy;
      }
    }
  }
  if (minmu == FLT_MAX) minmu = 0.;
  *length = (float)per;
  return(minmu);
}

float dist2(p,q)
float *p,*q;
{
  return((float) hypot((double)(*p-*q),(double)(*(p+1)-*(q+1))));
}


/* allocate memory for boundaries data */
void pixels_and_data(tree,NormofDu,image,prec,tabsaddles,sumsqper)
Shapes tree;
Fimage NormofDu,image;
float **tabsaddles,*sumsqper;
int *prec;
{
  Shape s;
  int i,k;
  Mydata data;
  Flist l;

  /* fill the field pixels of each shape*/
  mwdebug("Fill pixels field of each shape...");
  flst_pixels(tree);
  mwdebug("done.\n");
  /* compute the boundaries and allocate data field */
  mwdebug("Boundary length computation ...");
  l = mw_new_flist();
  *sumsqper = 0.;
  for(i=0;i<tree->nb_shapes;i++){
    s = tree->the_shapes+i;
    s->data = malloc(sizeof(struct mydata));
    if(!s->data) mwerror(FATAL,1,"Cannot allocate memory.\n");
    data = (Mydata) s->data;
    data->nfa = FLT_MAX;
    data->visit = 0;
    if(i>0) {
      flstb_boundary(prec,image,tree,s,NULL,l,tabsaddles);
      data->min_contrast = min_contrast(l,&data->length,NormofDu);
      *sumsqper += data->length*data->length;
    } else data->min_contrast = 0.;
    data->p_contrast = 0.;
    data->meaningful = (char) 0;
    data->notmax = (char) 0;
    data->max = (char) 0;
    data->type = (char) 0;
    data->ndetect = 0;
    data->histo = NULL;
    if(s->parent) s->removed = (char) 1;
  }
  mwdebug("done.\n");
  mw_delete_flist(l);
  fflush(stdout);
}


float image_max(image)
Fimage image;
{
  int i;
  float output;
  
  output = image->gray[0];
  for(i=0;i<image->ncol*image->nrow;i++)
    output = MAX(output,image->gray[i]);
  return(output);
}


/* compute gradient histogram  in a shape */
Fsignal shape_histo(shape,NormofDu,size,scale)
Shape shape;
Fimage NormofDu;
int size;
float scale;
{
  Point_plane p;
  int i,ncol,nrow,k,nbpt=0,ox,oy,x,y;
  float value,*tabvalues;
  
  Fsignal histo;
  
  ncol = NormofDu->ncol;
  nrow = NormofDu->nrow;
  
  histo = mw_change_fsignal(NULL,size);
  if(!histo) mwerror(FATAL,1,"Not enough memory.\n");
  mw_clear_fsignal(histo,0.);
  histo->scale = scale;
  /* compute histogram of gradient in shape */
  for(i=0,p=shape->pixels;i<shape->area;i++,p++){
    value = NormofDu->gray[p->x+ncol*p->y]/histo->scale;
    k = (int)value;
    if(k>=0) histo->values[k]+=1.;
  }
  
  /* because Du != 0 on level lines  */
  histo->values[0] = 0.;     
  return(histo);
}

/* compute log10 of repartition function from a histogram */
void update_repart(local_histo,local_repart)
Fsignal local_histo,local_repart;
{
  int nbins,i;

  nbins = local_histo->size;
  /* compute log10 of Proba(|Du|>x) */
  local_repart->values[nbins-1] = local_histo->values[nbins-1]+1;
  for(i=nbins-1;i>0;i--)
    local_repart->values[i-1] = local_repart->values[i]
      +local_histo->values[i-1];
  for(i=nbins-1;i>=0;i--) 
    local_repart->values[i] /= local_repart->values[0];
/*   affiche_repart(local_repart);  */
  for(i=0;i<nbins;i++)
    local_repart->values[i] = (float) log10(local_repart->values[i]);
}



float logH(mu,local_repart)
float mu;
Fsignal local_repart;
{
  int i;
  
  i = (int) (mu/local_repart->scale);
  if(i>= local_repart->size-1) i = local_repart->size-1;
  return(local_repart->values[i]);
}


/* compute IN PLACE the difference of two signals and take the positive part*/
void signal_pos_diff(sig1,sig2)
Fsignal sig1,sig2;
{
  int i;
  
  for(i=0;i<sig1->size;i++){
    sig1->values[i] -= sig2->values[i];
    if(sig1->values[i]<0) sig1->values[i] = 0.;
  }
}

void update_root_histo(local_root,local_histo,NormofDu)
Shape local_root;
Fsignal local_histo;
Fimage NormofDu;
{
  Shape s;
  Mydata sdata;
  int ncol,nrow,x,y,ox,oy,i,k;
  float *tabvalues,value;
  
  ncol = NormofDu->ncol;
  nrow = NormofDu->nrow;
  
  s = mw_get_first_child_shape(local_root);
  while(s){
    sdata = (Mydata) s->data;
    signal_pos_diff(local_histo,sdata->histo);
    s = mw_get_next_sibling_shape(s);
  }
  
  /*remove root boundary points from histogram (unless shape is absolut root)*/
  if(local_root->boundary){
    tabvalues = local_root->boundary->values;
    ox = MIN(ncol-1,MAX(0,(int) rint(tabvalues[0])-1));  
    oy = MIN(nrow-1,MAX(0,(int) rint(tabvalues[1])-1));
    for(i=0;i<local_root->boundary->size;i++){
      x = MIN(ncol-1,MAX(0,(int) rint(tabvalues[2*i])-1));  
      y = MIN(nrow-1,MAX(0,(int) rint(tabvalues[2*i+1])-1));
      if(x!= ox || y!=oy){
	ox = x; oy = y;
	value = NormofDu->gray[x+ncol*y]/local_histo->scale;
	k = (int)value;
	local_histo->values[k]-=1.;
	if(local_histo->values[k]<0.) local_histo->values[k]=0.;
      } 
    }
  }
}


/* compute nfa of a given shape boundary */
int update_mydata(s,type,local_root,local_repart,Global,new_open)
Shape s,local_root;
char type;
Fsignal local_repart;
struct myglobal Global;
int *new_open;
{ 
  Shape t;
  float mu,length,threshold;
  struct mydata *sdata, *tdata;
  int detect=0;
  
  if (!s) return;

  sdata = (Mydata) s->data;
  mu = sdata->min_contrast;
  threshold = Global.threshold;
  /* update nfa only if shape is not already as maximal or meaningful and not
     maximal.
     This is only useful for a local research where a curve may be 
     visited several times. 
  */
  if(!sdata->notmax && !sdata->max){     
    if (s != local_root) {
      
      /*** compute nfa (-log10) ***/
      sdata->p_contrast = logH(mu,local_repart);
      sdata->nfa = sdata->length*sdata->p_contrast/2.;
      sdata->visit++;
      
      /* update type */
      if (sdata->nfa < threshold) {
	sdata->type = s->inferior_type;
	sdata->meaningful = (char) 1;
	detect = 1;
	/* detect new open meaningful boundary; 
	   in the local research, we first deal with all open shapes before
	   detecting closed one. This is because the interior of an open shape 
	   is not clearly defined;
	   useful for local research only*/
	if(s->open) *new_open = 1;
      }
      else sdata->type = type;
    } else {
      sdata->nfa = threshold;
      sdata->type = s->inferior_type;
    }
  }
  
  /* visit children and update bestnfa_inf from children */
  sdata->ndetect = 0;
  if(s == local_root || !sdata->max)
    for (t=s->child;t;t=t->next_sibling){
      detect += update_mydata(t,sdata->type,local_root,local_repart,
			      Global,new_open);
      tdata = (struct mydata *)t->data;
      if (tdata->nfa<threshold) {
	/* useful for local research only:
	   increment the number of newly detected in s*/