ml_extract.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- Commande MegaWave -----------------------------*/
/* mwcommand
   name = {ml_extract};
   version = {"8.1"};
   author = {"Georges Koepfler"};
   function = {"Extract morpho_lines of image"};
   usage = {
   'L':level->level
       "for the value `level' compute morpho_lines (float)",
   'l':levels->levels
       "for each value in `levels' compute morpho_lines (Fsignal)",
   'o':[ml_opt=0]->opt [0,4]
       "select the content of the morpho_lines (min-max value)",
   'd':c_out<-c_out
       "draw the output mimage using ll_draw",
   'm'->m_flag    
       "optimize memory occupation during extraction",
   image_in->image_org
       "original image",
   mimage<-m_image
       "mimage with morpho_lines of image_in"
       };
*/

#include <stdio.h>
#include <math.h>

#define NDEBUG             /* comment this line out to enable assert() */
#include <assert.h>
#include "mw.h"

extern Cimage ml_draw();

#define INSIDE(A)   ((minvalue<=(A))&&((A)<=maxvalue))
#define OUTSIDE(A)  (!INSIDE(A))

void
produce_HV(im,NL,NC,minvalue,maxvalue,H,V)
float **im,minvalue,maxvalue;
unsigned int NL,NC;
unsigned char **H, **V;
{
  unsigned int l,c;

  for(l=1;l<NL;l++) for(c=1;c<NC;c++) {
    if(INSIDE(im[l][c])) {
      if(OUTSIDE(im[l][c-1]))  V[l][c-1]=1; 
      if(OUTSIDE(im[l-1][c]))  H[l-1][c]=1;
    }
    else {
      if(INSIDE(im[l][c-1])) V[l][c-1]=1; 
      if(INSIDE(im[l-1][c])) H[l-1][c]=1;
    }
  }
  /* oops: almost forgot these */
  for(l=1;l<NL;l++) {
    if(INSIDE(im[l][0])) 
      {if(OUTSIDE(im[l-1][0]))  H[l-1][0]=1;}
    else
      {if(INSIDE(im[l-1][0])) H[l-1][0]=1;}
  }
  for(c=1;c<NC;c++) {
    if(INSIDE(im[0][c]))
      {if(OUTSIDE(im[0][c-1]))  V[0][c-1]=1;}
    else
      {if(INSIDE(im[0][c-1])) V[0][c-1]=1;}
  }
}

unsigned long
count_X(H,V,NL,NC)
unsigned char **H,**V;
unsigned int NL,NC;
{
  unsigned int l,c,sum;
  unsigned long count_points=0L;

  for(l=0;l<NL-1;l++) for(c=0;c<NC-1;c++) {
    sum=V[l][c]+H[l][c]+V[l+1][c]+H[l][c+1];
    if(sum!=0) {
      assert((sum==2)||(sum==4));
      count_points++;
      if(sum==4) count_points++;
    }
  }
  for(c=0;c<NC-1;c++) {
    if(V[0][c]==1)    count_points++;
    if(V[NL-1][c]==1) count_points++;
  }
  for(l=0;l<NL-1;l++) {
    if(H[l][0]==1)    count_points++;
    if(H[l][NC-1]==1) count_points++;
  }
  return(count_points);
}

Morpho_line
produce_lline(old_lline,minvalue,maxvalue,current_point,current_type,open)
Morpho_line old_lline;
float minvalue,maxvalue;
Point_curve current_point;
Point_type  current_type;
unsigned char open;
{
  Morpho_line lline=NULL;

  if(!(lline=mw_change_morpho_line(lline)))
	  mwerror(FATAL,1,"Not enough memory for produce_lline.");
  lline->previous=NULL;
  lline->next=old_lline;
  lline->minvalue=minvalue;
  lline->maxvalue=maxvalue;
  lline->open=open;
  lline->first_point=current_point;
  lline->first_type=current_type;
  return(lline);
}


void
follow_open_line(NL,NC,H,V,ll,cc,sum,p,t)
unsigned int NL,NC;
unsigned char **H,**V;
int ll,cc,sum;
Point_curve *p;
Point_type  *t;
{
  Point_curve current_point=*p;
  Point_type  current_type=*t;

  while(sum!=0) {
    assert((sum==1)||(sum==3));
    current_point->x=cc+1;
    current_point->y=ll+1;
    current_point=current_point->next;
    current_type->type=0;                /* point inside image */
    current_type=current_type->next;
    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=((ll<0)||(cc<0)||(ll==NL-1)||(cc==NC-1)) ?
      0:(V[ll][cc]+H[ll][cc]+V[ll+1][cc]+H[ll][cc+1]);
  } /* end while */
  assert((ll<0)||(cc<0)||(ll==NL-1)||(cc==NC-1));
  current_point->y=(ll<0)? 0:((ll==NL-1)? NL:ll+1);
  current_point->x=(cc<0)? 0:((cc==NC-1)? NC:cc+1);
  *p=current_point->next;
  current_point->next=NULL;
  current_type->type=1;                    /* point on border of image */
  *t=current_type->next;
  current_type->next=NULL;
}

void
get_open_lines(im,NL,NC,minvalue,maxvalue,H,V,lline,p,t)
float **im, minvalue,maxvalue;
unsigned int NL,NC;
unsigned char **H,**V;
Morpho_line *lline;
Point_curve *p;
Point_type *t;
{
  Morpho_line old_lline;
  Point_curve current_point=*p;
  Point_type  current_type=*t;
  int sum,l,c;

  /* open lines starting in the first row */
  for(c=0;c<NC-1;c++) 
    if((V[0][c]==1)&&OUTSIDE(im[0][c])) { /* such that levelset at left */
      old_lline=*lline;
      *lline=produce_lline(old_lline,minvalue,maxvalue,current_point,current_type,1);
      if(old_lline!=NULL) old_lline->previous=*lline;
      current_point->x=c+1;
      current_point->y=0;
      current_point->previous=NULL;
      current_point=current_point->next;
      current_type->type=1;                /* point on border of image */
      current_type->previous=NULL;
      current_type=current_type->next;
      V[0][c]=0;
      sum=H[0][c]+V[1][c]+H[0][c+1];
      assert(sum!=0);
      follow_open_line(NL,NC,H,V,0,c,sum,&current_point,&current_type);
    }

  /* open lines starting in the last row */
  for(c=0;c<NC-1;c++) 
    if((V[NL-1][c]==1)&&INSIDE(im[NL-1][c])) { /* ...levelset at left */
      old_lline=*lline;
      *lline=produce_lline(old_lline,minvalue,maxvalue,current_point,current_type,1);
      if(old_lline!=NULL) old_lline->previous=*lline;
      current_point->x=c+1;
      current_point->y=NL;
      current_point->previous=NULL;
      current_point=current_point->next;      
      current_type->type=1;                /* point on border of image */
      current_type->previous=NULL;
      current_type=current_type->next;
      V[NL-1][c]=0;
      sum=H[NL-2][c+1]+V[NL-2][c]+H[NL-2][c];
      assert(sum!=0);
      follow_open_line(NL,NC,H,V,(int)NL-2,c,sum,&current_point,&current_type);
    }

  /* open lines starting in the first column */
  for(l=0;l<NL-1;l++)
    if((H[l][0]==1)&&INSIDE(im[l][0])) { /* ...levelset at left       */
      old_lline=*lline;
      *lline=produce_lline(old_lline,minvalue,maxvalue,current_point,current_type,1);
      if(old_lline!=NULL) old_lline->previous=*lline;
      current_point->x=0;
      current_point->y=l+1;
      current_point->previous=NULL;
      current_point=current_point->next;
      current_type->type=1;                /* point on border of image */
      current_type->previous=NULL;
      current_type=current_type->next;
      H[l][0]=0;
      sum=V[l+1][0]+H[l][1]+V[l][0];
      assert(sum!=0);
      follow_open_line(NL,NC,H,V,l,0,sum,&current_point,&current_type);
    }

  /* open lines starting in the last column */
  for(l=0;l<NL-1;l++)
    if((H[l][NC-1]==1)&&OUTSIDE(im[l][NC-1])) { /* ...levelset at left  */
      old_lline=*lline;
      *lline=produce_lline(old_lline,minvalue,maxvalue,current_point,current_type,1);
      if(old_lline!=NULL) old_lline->previous=*lline;
      current_point->x=NC;
      current_point->y=l+1;
      current_point->previous=NULL;
      current_point=current_point->next;
      current_type->type=1;                /* point on border of image */
      current_type->previous=NULL;
      current_type=current_type->next;
      H[l][NC-1]=0;
      sum=V[l][NC-2]+H[l][NC-2]+V[l+1][NC-2];
      assert(sum!=0);
      follow_open_line(NL,NC,H,V,l,(int)NC-2,sum,&current_point,&current_type);
    }
  *p=current_point;
  *t=current_type;
}

void
follow_closed_line(H,V,ll,cc,sum,p,t)
unsigned char **H,**V;
int ll,cc,sum;
Point_curve *p;
Point_type *t;
{
  Point_curve current_point=*p;
  Point_type  current_type=*t;

  while(sum!=0) {
    assert((sum==1)||(sum==3));