segct.c

image processing including fourier,wavelet,segmentation etc.

      || block_reg    ==NULL  ||  block_heap  ==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]   */
                                       /*            |         |              */
  image.energie.e= 0.0;                /*            S---------S              */
  image.energie.l= (long) image.PAS    /*      s[i,j+1]       s[i+1,j+1]      */
                      * image.nbbords; /*              b[2i,j]                */
                                       /*                                     */
                                       /***************************************/
  image.lambda=0L;
  image.heap=block_heap;
  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; 


  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 */
    RegCanalInit(image_org,i,j,regptr2->canal);
    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+=image.NC_REG;                     /* 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 */
      BorCanalInit(i,j,block_bor[l].canal);             /*       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   */
       BorCanalInit(i,j,block_bor[l].canal);                     /*   ^   */
       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*image.PAS;/*On place la frontiere dans la region droite*/
      block_som[l].y=j*image.PAS;/*           respectivement la region dessous*/
/* Not necessary in this version...
      SomCanalInit(i,j,block_som[l].canal);
*/
      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;
    }
  }

  printf("   ...\n");  /* here we initialize the heap stack */

  for (j=n=0;j<image.gy;j++)
  for (i=0;i<image.gx;i++) {
    if (j!=image.gy-1) {
      l=(long) 2*i  +lbbor*j;
      init_heap(block_bor+l,n++);
    }
    if (i!=image.gx-1) {
      l=(long) 2*i+1+lbbor*j;
      init_heap(block_bor+l,n++);
     }
  }
  printf(" finished !\n");
}

/* free all the 'callocated' memory blocks */
void
Image_free(imageptr)
MODELEPTR imageptr;
{
  free((void*)imageptr->heap);
  free((void*)imageptr->b_data);
  free((void*)imageptr->b_reg);
  free((void*)imageptr->b_bor);
  free((void*)imageptr->b_som);
  free((void*)imageptr->b_borcnxe);
  free((void*)imageptr->b_libor1);
  free((void*)imageptr->b_libor2);
  free((void*)imageptr->b_lipix); 
}


/* Merging tools */

void RegMerge(reg,regvois,bordcom)/* On agglutine les regions reg et regvois */
REGIONPTR reg,regvois;            /* de bord commun bordcom,dans l'espace    */
BORDPTR bordcom;                  /* memoire de reg                          */
{
  unsigned short i;
  BORDCONNEXEPTR cnxeptr;
  static void LiBordsUnion(),ElimBordeSom(),DegreSommet(),ElimLiReg(),UnionBordCnxe();


  for (i=0;i<image.NC_REG;i++)      reg->canal[i] += regvois->canal[i];

  if (image.regpixel0_0==regvois)  image.regpixel0_0=reg;

  ElimLiReg(regvois);

  LiBordsUnion(reg,regvois,bordcom);

  ElimBordeSom(bordcom);

  image.nbregions--;image.nbbords--;
  subtract_heap(bordcom);

  for(cnxeptr=bordcom->cnxe;cnxeptr!=NULL;cnxeptr=cnxeptr->suiv) {
        DegreSommet(cnxeptr->sa);
        if (cnxeptr->sb != cnxeptr->sa) DegreSommet(cnxeptr->sb);
  }
  UnionBordCnxe(reg);
}


void ElimLiReg(regvois)             /* Elimine la region regvois de la liste */
REGIONPTR regvois;                  /* des regions,exit si deux regions      */
{                                   /* Utiliser avant LiBordsUnion()         */

  if (  ((regvois->Rprec==NULL)&&(regvois->Rsuiv->Rsuiv==NULL))  || 
        ((regvois->Rsuiv==NULL)&&(regvois->Rprec->Rprec==NULL))     ) {
    printf("\nThere is just one region which stays, no boundaries!\n");
    exit(0);
  }
  if (regvois->Rprec==NULL) {
    image.tete=regvois->Rsuiv;
    image.tete->Rprec=NULL;
  }
  else 
    if (regvois->Rsuiv==NULL)
      regvois->Rprec->Rsuiv=NULL;
    else {
      regvois->Rprec->Rsuiv=regvois->Rsuiv;
      regvois->Rsuiv->Rprec=regvois->Rprec;
    }
 }


void LiBordsUnion(reg,regvois,bordcom)/*  Fait l'union,suivant l'ordre de    */
REGIONPTR reg,regvois;        /*  parcours de la frontiere de reg+regvois,de */
BORDPTR  bordcom;             /*  la liste de bords de reg et regvois.       */
{                             /* En eliminant l'element pointant sur bordcom */
  LI_BORDSPTR li,L1,L2;

  /* Marquer bordcom dans la liste des bords de reg */
  if (reg->Lbords->bord==bordcom)   L1=reg->Lbords ;
  else     {
            for(li=reg->Lbords;li->suiv->bord!=bordcom;li=li->suiv);
            L1=li->suiv;
           }

  /* Les bords de regvois doivent pointer vers reg,marquer bordcom */
  for (li=regvois->Lbords;li!=NULL;li=li->suiv)
    {
      if (li->bord==bordcom) L2=li;
      if (li->bord->rg==regvois)  li->bord->rg=reg ; 
      else                        li->bord->rd=reg ;
    }

  /* Mise en place ,dans l'ordre, des bords de reg (sens de rot. +).*/
  /* Si  reg->Lbords = a-L1-b  et regvois->Lbords = c-L2-d  ,       */
  /*  alors le resultat est   reg->Lbords = a-d-c-b .               */
  if (reg->Lbords==L1)
     if (regvois->Lbords==L2)
         if (L1->suiv==NULL)
             if (L2->suiv==NULL) {      /*   L1-0    et   L2-0   ,li=0        */
                li=NULL;
             }
             else                {      /*   L1-0    et   L2-d   ,li=d        */
                li=L2->suiv;
             }
         else
             if (L2->suiv==NULL) {      /*   L1-b    et   L2-0   ,li=b        */
                li=L1->suiv;
             }
             else                {      /*   L1-b    et   L2-d   ,li=d-b      */
                for (li=L2;li->suiv!=NULL;li=li->suiv)    ;
                li->suiv=L1->suiv;
                li=L2->suiv;
             }
     else
         if (L1->suiv==NULL)
             if (L2->suiv==NULL) {      /*   L1-0    et   c-L2-0 ,li=c        */
                for (li=regvois->Lbords;li->suiv!=L2;li=li->suiv) ;
                li->suiv=NULL;
                li=regvois->Lbords;
             }
             else                {      /*   L1-0    et   c-L2-d ,li=d-c      */
                for (li=L2;li->suiv!=L2;li=li->suiv)
                     if (li->suiv==NULL) li->suiv=regvois->Lbords;
                li->suiv=NULL;
                li=L2->suiv;
             }
         else
             if (L2->suiv==NULL) {      /*   L1-b    et   c-L2-0 ,li=c-b      */
                for (li=regvois->Lbords;li->suiv!=L2;li=li->suiv) ;
                li->suiv=L1->suiv    ;
                li=regvois->Lbords;
             }
             else                {      /*   L1-b    et   c-L2-d  ,li=d-c-b   */
                for (li=L2;li->suiv!=L2;li=li->suiv)
                     if (li->suiv==NULL) li->suiv=regvois->Lbords;
                li->suiv=L1->suiv;
                li=L2->suiv;
             }
  else
     if (regvois->Lbords==L2)
         if (L1->suiv==NULL)
             if (L2->suiv==NULL) {      /*   a-L1-0  et   L2-0    ,li=a       */
                for (li=reg->Lbords;li->suiv!=L1;li=li->suiv)   ;
                li->suiv=NULL;
                li=reg->Lbords;
             }
             else                {      /*   a-L1-0  et   L2-d    ,li=a-d     */
                for (li=reg->Lbords;li->suiv!=L1;li=li->suiv);
                li->suiv=L2->suiv;
                li=reg->Lbords;
             }
         else
             if (L2->suiv==NULL) {      /*   a-L1-b  et   L2-0    ,li=a-b     */
                for (li=reg->Lbords;li->suiv!=L1;li=li->suiv) ;
                li->suiv=L1->suiv;
                li=reg->Lbords;
             }
             else                {      /*   a-L1-b  et   L2-d    ,li=a-d-b   */
                for (li=reg->Lbords;li->suiv!=NULL;li=li->suiv)  
                                      if (li->suiv==L1) li->suiv=L2->suiv;
                li->suiv=L1->suiv;
                li=reg->Lbords;
             }
     else
         if (L1->suiv==NULL)
             if (L2->suiv==NULL) {      /*   a-L1-0  et   c-L2-0  ,li=a-c     */
                for (li=reg->Lbords;li->suiv!=L2;li=li->suiv) 
                         if (li->suiv==L1) li->suiv=regvois->Lbords ;
                li->suiv=NULL;
                li=reg->Lbords;
             }
             else                {      /*   a-L1-0  et   c-L2-d  ,li=a-d-c   */
                for (li=reg->Lbords;li->suiv!=L2;li=li->suiv) {
                     if (li->suiv==L1)   li->suiv=L2->suiv;
                     if (li->suiv==NULL) li->suiv=regvois->Lbords;
                }
                li->suiv=NULL;
                li=reg->Lbords;
             }
         else
             if (L2->suiv==NULL) {      /*   a-L1-b  et   c-L2-0  ,li=a-c-b   */
                for (li=reg->Lbords;li->suiv!=L2;li=li->suiv) 
                         if (li->suiv==L1) li->suiv=regvois->Lbords   ;
                li->suiv=L1->suiv;
                li=reg->Lbords;
             }