one_levelset.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- Commande MegaWave -----------------------------*/
/* mwcommand
   name = {one_levelset};
   version = {"1.01"};
   author = {"Georges Koepfler"};
   function = {"Get boundaries of level set, using a simplified merging criterion in the 'well-known' segmentation algorithm"};
   usage = {
   'l':[level=127.0]->level     [0.0,4e4]
       "pixels <=`level' (float) belong to the level set, default 127 ",
   'b':boundary<-cb    
       "output boundary of levelset, file cimage formated", 
   'p':polygons<-pb
       "output boundary of levelset, file fpolygons formated",
   'G':f_levelset<-fu
       "output levelset with gray-`level', file fimage formated",
   'B':b_levelset<-cu
       "output levelset b/w, file cimage formated",
   fimage->image_org
       "original image"
       };
*/
/*--- MegaWave2 - Copyright (C) 1994 Jacques Froment. All Rights Reserved. ---*/
#include <stdio.h>
#include <assert.h>
#include <math.h>
 
/* Include always the MegaWave2 Library */
#include "mw.h"
                                      /* Define the max,min of two quantities */
#define   MAX(A,B)          ((A) > (B)) ?  (A) : (B)
#define   MIN(A,B)          ((A) < (B)) ?  (A) : (B)

#define VRAI        1                     /*      Valeurs                     */
#define FAUX        0                     /*          logiques                */

/* #define NC_REG      1 !                   Nombre de canaux pour une region */
#define NC_BOR      1                     /* Nombre de canaux pour un bord    */
#define NC_SOM      0                     /* Nombre de canaux pour un sommet  */

/* Structure definitions for segmentation */

 typedef struct li_pixels {
    char              direction;
    unsigned short    longueur;
    struct li_pixels *suiv;
 } LI_PIXELS,*LI_PIXELSPTR;

 typedef struct li_bords {
     struct bord      *bord;
     struct li_bords  *suiv;
 } LI_BORDS,*LI_BORDSPTR;

 typedef struct region {
     float              *canal;    /* in this algorithm only the gray level */
     LI_BORDSPTR        Lbords;    /* is of interest !                      */
     struct region     *Rprec,*Rsuiv;
 } REGION,*REGIONPTR;

 typedef struct bord {                     /********************************/
     float              canal;             /*       sa          lab   lba  */
     REGIONPTR          rg,rd;             /*        |           |     ^   */
     struct bord_cnxe  *cnxe;              /*    rg  |  rd   ;   |     |   */
 } BORD,*BORDPTR;                          /*        |           |     |   */
                                           /*       sb           v     |   */
                                           /********************************/
 typedef struct bord_cnxe {
     struct sommet     *sa,*sb;
     LI_PIXELSPTR       lab,lba;
     struct bord_cnxe  *suiv;
 } BORDCONNEXE,*BORDCONNEXEPTR;

 typedef struct sommet {
/*Not used in this version
    float             canal[NC_SOM];
*/
    unsigned short    x,y;
    LI_BORDSPTR       Lbords;
 } SOMMET,*SOMMETPTR;

 typedef struct modele {
      unsigned short   gx,gy;            /* Dimensions de la grille           */
      unsigned long    nbregions;        /* Nombre de regions                 */
      unsigned long    nbbords,nbsommets;/* Nombre de bords,resp. de sommets  */
      REGIONPTR        regpixel0_0;      /* Region qui contient l'origine     */
      REGIONPTR        tete;             /* Tete de la liste image            */
      /* pointers to the 'callocated' memory blocks                           */
      float *          b_data;           /* float-region-canal-block          */
      REGIONPTR        b_reg;            /* region-block                      */
      BORDPTR          b_bor;            /* bords-block                       */
      SOMMETPTR        b_som;            /* sommets-block                     */
      BORDCONNEXEPTR   b_borcnxe;        /* connected comp-block              */
      LI_BORDSPTR      b_libor1,b_libor2;/* bords-liste-block                 */
      LI_PIXELSPTR     b_lipix;          /* pixels-liste-block                */
} MODELE,*MODELEPTR;


/* DECLARATION DE LA VARIABLE GLOBALE */

MODELE image;                               /* Image sous la forme du modele */

/* Channel initilization functions */

void InitPixel(pixelptr,dir)        /* Met un segment de longeur 1     et de */
LI_PIXELSPTR pixelptr;              /*  direction 'dir' dans  'pixelptr'.    */
char dir;
{
  pixelptr->direction=dir;
  pixelptr->longueur=1;
  pixelptr->suiv=NULL;
}

/* Initilization of the intern image model */

void Estimation()                          /* Estimations occupation memoire */
{
  unsigned long  simage,stotal,total,
                 sregion,sbord,sbordcnxe,ssommet,sli_bords,sli_pixels,
                 sdata,
                 nbregions,nbbords,nbsommets,systeme;

  sregion     =  sizeof(REGION);
  sbord       =  sizeof(BORD);
  sbordcnxe   =  sizeof(BORDCONNEXE);
  ssommet     =  sizeof(SOMMET);
  sli_bords   =  sizeof(LI_BORDS);
  sli_pixels  =  sizeof(LI_PIXELS);
  sdata       =  sizeof(float);

  nbregions =  (image.gx*image.gy);
  nbbords   =  (image.gx*(image.gy-1)+image.gy*(image.gx-1));
  nbsommets =  ((image.gx+1)*(image.gy+1)-4);

  simage=   nbregions  * (sregion+sdata)
         +  nbbords    * (sbord+sbordcnxe)
         +  nbsommets  *     ssommet
         +  2*nbbords  * (2*sli_bords+sli_pixels)   ;

  systeme= 200000.0;  /* On utilise approximativement 200 kO */
                      /* pour la reservation par 7 'calloc()'*/

  stotal=systeme+image.gx*image.gy*sizeof(float);
  total= simage + stotal;
  printf("\nEstimation for memory occupation: %d Mbytes\n",total>>20);
  return;
}
 
void 
Initialisation(image_org)       
Fimage  image_org;     /* On cree huit blocs pour reserver de la memoire      */
{                      /* aux regions,bords et sommets,... ; Grace a un codage*/
                       /* des places memoire on peut facilement initialiser   */
  unsigned short i,j;  /* les pointeurs des differentes structures.           */
  unsigned long  l,n,
                 dbreg,dbbor,dbsom,dblibor,dblipix,   /* Dimensions des blocs */
                 lbreg,lbbor,lbsom;                   /* Largeur des blocs    */

  REGIONPTR       block_reg,                /* Espace pour les regions        */
                  regptr1,regptr2;
  BORDPTR         block_bor;                /* Espace pour les bords          */
  SOMMETPTR       block_som;                /* Espace pour les sommets        */
  BORDCONNEXEPTR  block_borcnxe;            /* Espace pour les comp. connexes */
  LI_BORDSPTR     block_libor1,block_libor2,/* Espace pour el. liste  de bords*/
                  liborptr1,liborptr2;
  LI_PIXELSPTR    block_lipix;              /* Espace pour liste de pixels    */
  float *         block_data;               /* Espace pour les canaux         */

  static void InitPixel();


  printf("\nInitialization\n");
  lbreg=image.gx;     dbreg=lbreg*image.gy;
  lbbor=2*image.gx-1; dbbor=lbbor*image.gy;
  lbsom=image.gx+1;   dbsom=lbsom*(image.gy+1);
  dblibor=dblipix=2*((image.gx-1)*image.gy+image.gx*(image.gy-1))+1L;

  block_reg     = (REGIONPTR     ) calloc((long)  dbreg  ,sizeof(REGION     ));
  block_data    = (float   *     ) calloc((long)  dbreg  ,sizeof(float      ));
  block_bor     = (BORDPTR       ) calloc((long)  dbbor  ,sizeof(BORD       ));
  block_som     = (SOMMETPTR     ) calloc((long)  dbsom  ,sizeof(SOMMET     ));
  block_borcnxe = (BORDCONNEXEPTR) calloc((long)  dbbor  ,sizeof(BORDCONNEXE));
  block_libor1  = (LI_BORDSPTR   ) calloc((long)  dblibor,sizeof(LI_BORDS   ));
  block_libor2  = (LI_BORDSPTR   ) calloc((long)  dblibor,sizeof(LI_BORDS   ));
  block_lipix   = (LI_PIXELSPTR  ) calloc((long)  dblipix,sizeof(LI_PIXELS  ));

  if (   block_lipix  ==NULL  ||  block_libor2==NULL  || block_libor1==NULL
      || block_borcnxe==NULL  ||  block_som   ==NULL  || block_bor   ==NULL
      || block_reg    ==NULL  ||  block_data  ==NULL ) 
     mwerror(FATAL,1,"\n Not enough memory for initialisation!!\n");

                                       /***************************************/
                                       /*  CODAGE DES REGIONS,BORDS,SOMMETS   */
  image.tete=block_reg;                /*     DANS LES BLOCS RESPECTIFS       */
  image.regpixel0_0=image.tete;        /*                                     */
  image.nbregions=image.gx*image.gy;   /*             b[2i,j-1]               */
  image.nbbords=image.gx*(image.gy-1)+ /*        s[i,j]       s[i+1,j]        */
                (image.gx-1)*image.gy; /*            S---------S              */
  image.nbsommets=(image.gx+1)         /*            |         |              */
                    *(image.gy+1)-4;   /*  b[2i-1,j] |  r[i,j] |  b[2i+1,j]   */
                                       /*            |         |              */
                                       /*            S---------S              */
                                       /*      s[i,j+1]       s[i+1,j+1]      */
                                       /*              b[2i,j]                */
                                       /*                                     */
                                       /***************************************/

  image.b_data=block_data;
  image.b_reg=block_reg;
  image.b_bor=block_bor;
  image.b_som=block_som;
  image.b_borcnxe=block_borcnxe;
  image.b_libor1=block_libor1;
  image.b_libor2=block_libor2;
  image.b_lipix=block_lipix;

  for(l=0;l<image.gx*image.gy;l++)     /* copy the gray-level */
   block_data[l]=image_org->gray[l];

  printf("   ...of the regions,   initial number of regions    %ld\n",image.nbregions);  
  regptr1=image.tete;  
  for (j=0;j<image.gy;j++)
  for (i=0;i<image.gx;i++) {
    regptr1->Rprec = ((i==0)&&(j==0))    ?    NULL    :   regptr2  ;
    regptr2=regptr1;
    liborptr1=liborptr2=NULL;
    regptr2->canal=block_data;           /* avoid computations */
    if (i!=image.gx-1)  {                /* Bord droit */
      liborptr2= block_libor1 ++ ;
      liborptr2->bord=( block_bor + ((long) (2*i+1) + lbbor* j) );
      liborptr1=liborptr2;
    }
    if (j!=0)  {                         /* Bord dessus */
      if (liborptr2==NULL)     liborptr2= block_libor1 ++ ;
      else   {liborptr2->suiv= block_libor1 ++ ;liborptr2=liborptr2->suiv;}
      liborptr2->bord=( block_bor + ((long) ( 2*i ) + lbbor*(j-1)) );
      if (liborptr1==NULL) liborptr1=liborptr2;
    }
    if (i!=0)  {                         /* Bord gauche */
      if (liborptr2==NULL)     liborptr2= block_libor1 ++ ;
      else   {liborptr2->suiv= block_libor1 ++ ;liborptr2=liborptr2->suiv;}
      liborptr2->bord=( block_bor + ((long) (2*i-1) + lbbor* j) );
      if (liborptr1==NULL) liborptr1=liborptr2;
    }
    if (j!=image.gy-1)  {                /* Bord dessous */
      if (liborptr2==NULL)     liborptr2= block_libor1 ++ ;
      else   {liborptr2->suiv= block_libor1 ++ ;liborptr2=liborptr2->suiv;}
      liborptr2->bord=( block_bor + ((long) ( 2*i ) + lbbor* j) );
      if (liborptr1==NULL) liborptr1=liborptr2;
    }
    liborptr2->suiv=NULL;
    regptr2->Lbords=liborptr1;
    regptr1=  ((i==image.gx-1) && (j==image.gy-1))     ?
                          NULL      :    (block_reg + ((long) j*lbreg+i+1));
    regptr2->Rsuiv=regptr1;
    block_data++;                               /* update the data pointer */
  }


  printf("   ...of the boundary,  initial number of frontiers  %ld\n",image.nbbords);
  for (j=0;j<image.gy;j++)
  for (i=0;i<image.gx;i++) {                            /*******************/
    if (j!=image.gy-1) {                                /*       rg        */
      l=(long) 2*i  +lbbor*j;                           /* sb --------> sa */
      block_bor[l].canal=1.0;                           /*       rd        */
      block_bor[l].rg=block_reg + (long) i + lbreg* j ; /*******************/
      block_bor[l].rd=block_reg + (long) i + lbreg*(j+1); 
      block_bor[l].cnxe= block_borcnxe + l;                    /**************/
      block_bor[l].cnxe->sa=block_som +(long)(i+1)+lbsom*(j+1);/*     'h'    */ 
      block_bor[l].cnxe->sb=block_som +(long)  i  +lbsom*(j+1);/* 'g' SOM 'd'*/
      block_bor[l].cnxe->lab= block_lipix ++;                  /*     'b'    */
      InitPixel(block_bor[l].cnxe->lab,'g');                   /**************/
      block_bor[l].cnxe->lba= block_lipix ++;
      InitPixel(block_bor[l].cnxe->lba,'d');
      block_bor[l].cnxe->suiv=NULL;
    }
    if (i!=image.gx-1) {                                         /*********/
       l=(long) 2*i+1+lbbor*j;                                   /*  sa   */
       block_bor[l].canal=1.0;                                   /*   ^   */
       block_bor[l].rg=block_reg + (long) i   +lbreg* j ;        /*   |   */
       block_bor[l].rd=block_reg + (long)(i+1)+lbreg* j ;        /* rg|   */
       block_bor[l].cnxe= block_borcnxe + l;                     /*   |rd */
       block_bor[l].cnxe->sa=block_som + (long)(i+1)+lbsom* j   ;/*   |   */
       block_bor[l].cnxe->sb=block_som + (long)(i+1)+lbsom*(j+1);/*  sb   */
       block_bor[l].cnxe->lab= block_lipix ++;                   /*********/
       InitPixel(block_bor[l].cnxe->lab,'b');
       block_bor[l].cnxe->lba= block_lipix ++;
       InitPixel(block_bor[l].cnxe->lba,'h');
       block_bor[l].cnxe->suiv=NULL;
    }
  }

  printf("   ...of tips,          initial number of tips       %ld\n",image.nbsommets);
  for (j=0;j<=image.gy;j++)
  for (i=0;i<=image.gx;i++) {
    liborptr1=liborptr2=NULL;
    l=(long) i+lbsom*j;
    if(((i!=0)||(j!=0&&j!=image.gy))&&((i!=image.gx)||(j!=0&&j!=image.gy))) {
      block_som[l].x=i;          /*On place la frontiere dans la region droite*/
      block_som[l].y=j;          /*            respectivement la region desous*/
      if ((i!= 0)&&(i!=image.gx)&&(j!=image.gy))  {      /* Bord en dessous */
        liborptr2= block_libor2 ++;
        liborptr2->bord=block_bor + (long)(2*i-1) + lbbor* j ;
        liborptr1=liborptr2;
      }
      if ((i!=image.gx)&&(j!= 0)&&(j!=image.gy))  {      /* Bord a droite   */
        if (liborptr2==NULL)     liborptr2= block_libor2 ++;
        else   {liborptr2->suiv= block_libor2 ++;liborptr2=liborptr2->suiv;}
        liborptr2->bord=block_bor + (long)( 2*i ) + lbbor*(j-1);
        if (liborptr1==NULL) liborptr1=liborptr2;
      }
      if ((i!= 0)&&(i!=image.gx)&&(j!= 0))  {      /* Bord en dessus  */
        if (liborptr2==NULL)     liborptr2= block_libor2 ++;
        else   {liborptr2->suiv= block_libor2 ++;liborptr2=liborptr2->suiv;}
        liborptr2->bord=block_bor + (long)(2*i-1) + lbbor*(j-1);
        if (liborptr1==NULL) liborptr1=liborptr2;
      }
      if ((i!= 0)&&(j!= 0)&&(j!=image.gy))  {      /* Bord a gauche   */
        if (liborptr2==NULL)     liborptr2= block_libor2 ++;
        else   {liborptr2->suiv= block_libor2 ++;liborptr2=liborptr2->suiv;}
        liborptr2->bord=block_bor + (long) 2*(i-1) + lbbor*(j-1);
        if (liborptr1==NULL) liborptr1=liborptr2;
      }
      liborptr2->suiv=NULL;
      block_som[l].Lbords=liborptr1;