segct.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- Commande MegaWave -----------------------------*/
/* mwcommand
   name = {segct};
   version = {"1.21"};
   author = {"Georges Koepfler"};
   function = {"Region-Growing method using the energy model of Mumford and Shah with piecewise constant approximation function"};
   usage = {
   'S':[size_of_grid=1]->sgrid
       "size of initilization grid (int), default 1 ",
   'N':nb_of_regions->nb_of_regions
       "number of desired regions (int) ",
   'L':lambda->lambda
       "value of final scale parameter (float) of last 2-normal segmentation",
   'c':curves<-curves
       "output boundary set in curves format",
   'r':reconstruction<-u
       "output piecewise constant reconstruction",
   f_nb_of_regions<-f_nb_of_regions
       "final number of regions",
   f_lambda<-f_lambda
       "final lambda value",
   fimage->image_org
       "original image",
   boundary<-segct
       "b/w image of boundary set "
          };
*/
/*--- MegaWave2 - Copyright (C) 1994 Jacques Froment. All Rights Reserved. ---*/
#include <stdio.h>
 
/* Include always the MegaWave2 Library */
#include "mw.h"


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

/* #define NC_REG      2                     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  */

/* Structures Definition */

 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 heap {
     float              lambda;
     struct bord       *bord;
 } HEAP,*HEAPTR; 

 typedef struct region {
     float             *canal;
     LI_BORDSPTR        Lbords;
     struct region     *Rprec,*Rsuiv;
 } REGION,*REGIONPTR;


 typedef struct bord {                    
     unsigned long      heap;              /********************************/
     float              canal[NC_BOR];     /*       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 energie {
     long  l;           /* Longeur des bords                                  */
     float e;           /* g etant l'image, l'energie elastique = S( (u-g)^2 )*/
 } ENERGIE;


 typedef struct modele {
      int              PAS;              /* value of sgrid                    */
      unsigned short   NC_REG;           /* nombre de canaux                  */
      float            lambda;           /* Echelle de la seg. courante       */
      unsigned short   gx,gy;            /* Dimensions de la grille           */
      unsigned long    nbregions;        /* Nombre de regions                 */
      unsigned long    nbbords,nbsommets;/* Nombre de bords,resp. de sommets  */
      ENERGIE          energie;          /* Energie elastique et long. l(B)   */
      REGIONPTR        regpixel0_0;      /* Region qui contient l'origine     */
      REGIONPTR        tete;             /* Tete de la liste image            */
      /* pointers to the 'callocated' memory blocks                           */
      HEAPTR           heap;             /* Base de la table ordonne des bords*/
      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;


/*DECLARATIONS DES VARIABLES GLOBALES*/

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


/* Channel initilization functions */


void RegCanalInit(image_org,i,j,canal)
                               /* Initialise les canaux associes aux regions */
Fimage image_org;
unsigned short i,j;
float *canal;
 {
  static float SomGris();

  canal[0]= (float) image.PAS*image.PAS;
  canal[1]= SomGris(image_org,i,j);
 }


float SomGris(image_org,i,j)   /* Calcule la somme de gris dans le rectangle */
Fimage image_org;
unsigned short i,j;            /*  (i*PAS,j*PAS)-( (i+1)*PAS-1,(j+1)*PAS-1 ).*/
{                              /* Determine l'energie elastique initale !!   */
  unsigned char k,l,c;
  float somme=0.0,sommeSpas;

  for (k=0;k<image.PAS;k++)
    for (l=0;l<image.PAS;l++) 
      somme+=(float)image_org->gray[(j*image.PAS+l)*image_org->ncol+i*image.PAS+k];
  return (somme);
}


void BorCanalInit(i,j,canal)     /* Initialise les canaux associes aux bords */
unsigned short i,j;              /* de la region [i,j]                       */
float canal[];                   /* Si i est pair c'est le bord du dessous,  */
{                                /* si i est impair c'est le bord de droite  */
  canal[0] = (float) image.PAS;
}


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


/* The macros and functions which manipulate the heapstack */
/* For  references see "Numerical recepies in C"           */
/* although this implementation is 'nicer'                 */

#ifndef u_long
#define u_long unsigned long
#endif
#define h_root           ((u_long)0)
#define h_add(A)         ((image.heap+(A)))
#define h_bord(A)        (image.heap[(u_long)(A)].bord)      /* type BORDPTR */
#define h_lambda(A)      (image.heap[(u_long)(A)].lambda)    /* type float   */
#define h_up(A)          ((u_long)((u_long)(A)-1)>>1)
#define h_left(A)        ((u_long)((u_long)(A)<<1)+1)
#define h_right(A)       ((u_long)((u_long)(A)+1)<<1)
#define bord_heap(B)     (((BORDPTR)(B))->heap)

void
init_heap(bord,n)
BORDPTR bord;
u_long n;
{
  u_long m;
  float lam;
  static float eval_lambda();
  
  lam=eval_lambda(bord);
  while(n>0 && (lam<h_lambda(m=h_up(n)))) {
    bord_heap(h_bord(m))=n;
    *h_add(n)=*h_add(m);
    n=m;
  }
  bord->heap=n;
  h_bord(n)=bord;
  h_lambda(n)=lam;
  return;
  }


void
update_heap(bord)
BORDPTR bord;
{
  u_long n,m0,m1;
  float lam;
  static float eval_lambda();

  n=bord_heap(bord);
  lam=eval_lambda(bord);
  if( n>0 && (lam<h_lambda(m0=h_up(n))))
    do{  /* bord is too light... it will bubble 'up' to the root */
      bord_heap(h_bord(m0))=n;
      *h_add(n)=*h_add(m0);
      n=m0;
    }while( n>0 && (lam<h_lambda(m0=h_up(n))) );
  else { /* bord is too heavy... it will sink to the ground */
    while( (m0=h_left(n))<image.nbbords) {
      if( ((m1=h_right(n))<image.nbbords) && (h_lambda(m1)<h_lambda(m0)))  m0=m1;
      if(h_lambda(m0)>lam) break;
      bord_heap(h_bord(m0))=n;
      *h_add(n)=*h_add(m0);
      n=m0;
    }
  }
  bord->heap=n;
  h_bord(n)=bord;
  h_lambda(n)=lam;
  return;
}

void
subtract_heap(bord)
BORDPTR bord;
{
  BORDPTR lbord;
  u_long n;

  n=bord_heap(bord);
  if(n<image.nbbords) {
    lbord=h_bord(image.nbbords); /* nbbords IS already decremented */
    bord_heap(lbord)=n;
    *h_add(n)=*h_add(image.nbbords);
    update_heap(lbord);
  }
  return;
}



/* Initilization of the intern image model */

void Estimation()                          /* Estimations occupation memoire */
{
  unsigned long  simage,stotal,total,
                 sregion,sbord,sbordcnxe,ssommet,sli_bords,sli_pixels,
                 sheap,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);
  sheap       =  sizeof(HEAP); /* no mechoui around ? */
  sdata       =  sizeof(float)*image.NC_REG;

  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+sheap) + 
              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.PAS*image.gy*image.PAS*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 */
                 dbheap,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    */
  HEAPTR          block_heap;               /* Espace pour heapstack          */
  float *         block_data;               /* Espace pour les canaux         */

  static void InitPixel(),RegCanalInit(),BorCanalInit();


  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;
  dbheap=(image.gx-1)*image.gy+image.gx*(image.gy-1);

  block_heap    = (HEAPTR        ) calloc((long)  dbheap ,sizeof(HEAP       ));
  block_reg     = (REGIONPTR     ) calloc((long)  dbreg  ,sizeof(REGION     ));
  block_data    = (float   *     ) calloc((long)  dbreg  ,
                                             image.NC_REG*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