infsup.c

image processing including fourier,wavelet,segmentation etc.

/* mwcommand 
name = {infsup};
version={"2.1"};
author = {"Frederic Guichard, Denis Pasquignon, Jacques Froment"};
function = {"InfSup scheme or median-median filtering"};
usage = {
'n':[Niter=2]->Niter[0,1000]      "number of iterations (default 2)",
'i':[deginf=0.0]->deginf[0,1]  "second level med (0..1, default: sup=0.)",
's':[degsup=1.0]->degsup[0,1]  "first level med (0..1, default: inf=1.)",
'a'->average "swap first and second level med and take the average",
input-> image  "input cimage",
fmovie -> fmovie "masks sequence",
output<- output   "filtered image"};
*/
/*----------------------------------------------------------------------
 v2.1: return void (L.Moisan)
----------------------------------------------------------------------*/

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

#include "mw.h"


int diff(a,b)
     unsigned char *a;
     unsigned char *b;
{
  if ((*a)<(*b)) {
    return (-1);
  }
  else {
    return (1);
  }
}

unsigned char TheInf(image,maskn,x,y)
     Cimage image;
     Fimage maskn;
     int x,y;
{
  int xd,yd,xf,yf,ntx,nty,Mdx,Mdy;
  int dx,dy;      /* taille de l'image */
  register int i,j;
  int aa,mm;      /* intermediaires de calcul */
  int l;
  unsigned char value;

  value=255;
  dx= image->ncol;
  dy= image->nrow;

  Mdx = maskn->ncol;
  Mdy = maskn->nrow;
  ntx= ((Mdx-1)/2);
  nty= ((Mdy-1)/2);

  xd=x-ntx;
  yd=y-ntx;
  xf=x+nty;
  yf=y+nty;
  if (xd<0) xd=0; else if (xf>=dx) xf=dx-1;
  if (yd<0) yd=0; else if (yf>=dy) yf=dy-1;
  aa= yd*dx;
  mm = -x+ntx+Mdx*(yd-y+nty);
  for(j=yd;j<=yf;j++) {
    for(i=xd;i<=xf;i++) {
      if (maskn->gray[mm+i]>0){ 
	if (value> (image-> gray[aa+i]))
	  value=image-> gray[aa+i];}
    }
    aa+=dx;
    mm+=Mdx;
  }
  
  return(value);
}

unsigned char TheSup(image,maskn,x,y)
     Cimage image;
     Fimage maskn;
     int x,y;
{
  int xd,yd,xf,yf,ntx,nty,Mdx,Mdy;
  int dx,dy;      /* taille de l'image */
  register int i,j;
  int aa,mm;      /* intermediaires de calcul */
  int l;
  unsigned char value;

  value=0;
  dx= image->ncol;
  dy= image->nrow;

  Mdx = maskn->ncol;
  Mdy = maskn->nrow;
  ntx= ((Mdx-1)/2);
  nty= ((Mdy-1)/2);

  xd=x-ntx;
  yd=y-ntx;
  xf=x+nty;
  yf=y+nty;
  if (xd<0) xd=0; else if (xf>=dx) xf=dx-1;
  if (yd<0) yd=0; else if (yf>=dy) yf=dy-1;
  aa= yd*dx;
  mm = -x+ntx+Mdx*(yd-y+nty);
  for(j=yd;j<=yf;j++) {
    for(i=xd;i<=xf;i++) {
      if (maskn->gray[mm+i]>0) { 
	if (value< (image-> gray[aa+i]))
	  value=image-> gray[aa+i];}
    }
    aa+=dx;
    mm+=Mdx;
  }
  
  return(value);
}



float  Histogram( image,maskn,x,y,Hist)
     Cimage image;
     Fimage maskn;
     int x,y;
     float *Hist;
{
  int xd,yd,xf,yf,ntx,nty,Mdx,Mdy;
  int dx,dy;      /* taille de l'image */
  register int i,j;
  int aa,mm;      /* intermediaires de calcul */
  int l;
  float sum;      /* sum of histogram  */

  dx= image->ncol;
  dy= image->nrow;

  Mdx = maskn->ncol;
  Mdy = maskn->nrow;
  ntx= ((Mdx-1)/2);
  nty= ((Mdy-1)/2);

  xd=x-ntx;
  yd=y-ntx;
  xf=x+nty;
  yf=y+nty;
  if (xd<0) xd=0; else if (xf>=dx) xf=dx-1;
  if (yd<0) yd=0; else if (yf>=dy) yf=dy-1;
  
  /* --- Histogramme est remis a zero */
  for(i=0;i<256;i++) Hist[i]=0.0;
  
  /* --- Puis on le calcule   */
  
  sum=0;
  aa= yd*dx;
  mm = -x+ntx+Mdx*(yd-y+nty);
  for(j=yd;j<=yf;j++) {
    for(i=xd;i<=xf;i++) {
      l=image-> gray[aa+i];
      Hist[l] += maskn->gray[mm+i];
      sum += maskn->gray[mm+i];
    }
    aa+=dx;
    mm+=Mdx;
  }
  
  return(sum);
}



/* Procedure InfSup  */

/* t1, t2 are the level of the median */
/* t1=1, t2=0, process a InfSup... */

ResolutionInfSup(image,A,G,t1,t2,movie,Nmask)
     Cimage image;
     Cimage A;              /* Image Buffer */
     unsigned char *G;      /* contains the values for each mask */     
     Fmovie movie;
     int Nmask;             /* Number of masks   */
     float t1,t2;
{
  Fimage maskn;
  int kx,ky;
  int i,j,l;
  register int k;
  int dx,dy;                /* size of the image  */
  int n,Num;
  int Mdx,Mdy;              /* size of a mask     */

  float Hist[257];          /* Histogramm         */
  float v, vmax;
  float presquezero;

  presquezero=0.00001;

  /* Init                                         */

  dx=image->ncol;
  dy=image->nrow;
  maskn=movie->first;
  Mdx=maskn->ncol;
  Mdy=maskn->nrow;
  
  if ((A == NULL) || (A->gray == NULL))
    mwerror(INTERNAL,1,"[ResolutionInfSup] Null buffer image !\n");

  if (G == NULL)
    mwerror(INTERNAL,1,"[ResolutionInfSup] Null G buffer !\n");

  Num= (int) ( (float) Nmask-1) * t2;

  /* For each point of the image : */

  for(j=0;j<dy;j++) {
    if ((j/10)*10==j) mwdebug("\t[%3.0f %%]\n",100.0 * (float) j/dy);
    for(i=0;i<dx;i++){

      
      maskn=movie->first;
      for(n=0;n<Nmask;n++) {                 /* For each mask */
      
	
	if ((t1>0.02)&&(t1<0.98))  /* -> median of level t1 */ 
	  {
	    vmax=Histogram( image,maskn,i,j,Hist);  /* Histogramm */
	    Hist[257]=100000.0;
	    vmax=vmax*t1;
	    v=presquezero; k=0; l=0;
	    while(Hist[k]<presquezero) k++; 
	    if (v<=vmax) {
	      while(v<=vmax) {   
		if (Hist[k]>presquezero) {
		  v+=Hist[k];
		}
		k++;
	      }
	      k--;
	      if (k==256) { k=255; while(Hist[k]<presquezero) k--;}
	    }
	  }
      
     

	/* if t1<0.02 -> inf */

	if (t1<=0.02) {
	  k=TheInf(image,maskn,i,j);
	}
	
	/* if t1>0.98 -> sup */
	
	if (t1>=0.98)
	  {
	    k=TheSup(image,maskn,i,j);
	  }
	
        /* we put the median into G */
	
	G[n]=k;
	
	maskn= maskn->next;
      
      }

      /* we estimate the median of level t2, on the */
      /* values of G                                */
      

      /* if (t2>0.02)&&(t2<0.98) -> median */

      if ((t2>0.02)&&(t2<0.98)) {
	qsort((void *) G, Nmask, sizeof(unsigned char), diff);
	k=G[Num];
      }
      else
	{
	  /* if t2<0.02 -> inf  */

	  if (t2<=0.02)
	    {
	      k=255;
	      for(l=0;l<Nmask;l++)
		if (k>G[l]) k=G[l];
	    }

	  /* if t2>0.98 -> sup */

	  if (t2>=0.98)
	    {
	      k=0;
	      for(l=0;l<Nmask;l++)
		if (k<G[l]) k=G[l];
	    }
	}
      /* We copy the result into the image A */

      A->gray[i+j*dx]=k;
    }
  }
  
  /* Finish */

  /* we copy the result in the original image */

  mw_copy_cimage(A,image);
}

/* ----------------------------------------------------------*/

/* The main procedure */

/* ----------------------------------------------------------*/


Cimage infsup(Niter,deginf,degsup,average,image,fmovie,output)
     int *Niter;
     float *deginf;
     float *degsup;
     char *average;
     Cimage image;
     Fmovie fmovie;
     Cimage output;
{
  int n,j,size,Nmask;
  Fimage i;
  float t1,t2;
  unsigned char *G=NULL;
  Cimage image_work=NULL;
  Cimage image_save=NULL;
  unsigned char *ptr1,*ptr2;

  /*------------------------  Number of images in the fmovie...  ------*/
  
  if ((fmovie -> first)==NULL) mwerror(FATAL,1,"empty masks movie !\n");
  Nmask=1;
  i = fmovie->first;
  while (i->next != NULL) { i=i->next; Nmask++; }
  mwdebug("Number of masks : %i \n",Nmask);

  /* ------------------------ Memory allocation for the buffers -----------*/

  G= (unsigned char *) calloc(Nmask+2 , sizeof(unsigned char));
  if (G == NULL) mwerror(FATAL, 1, "Not enough memory\n");

  image_work = mw_change_cimage(NULL,image->nrow,image->ncol);
  if (image_work == NULL) mwerror(FATAL, 1, "Not enough memory\n");  

  if ((output=mw_change_cimage(output,image->nrow,image->ncol)) == NULL)
    mwerror(FATAL, 1, "Not enough memory\n");  

  if (average)
    {
      image_save = mw_change_cimage(NULL,image->nrow,image->ncol);
      if (image_save == NULL) mwerror(FATAL, 1, "Not enough memory\n");
      mw_copy_cimage(image,image_save);
    }

  /* ----------------------- Main loop  ----------------------------------- */
  
  t1=(*degsup);
  t2=(*deginf);
  for(n=0;n<(*Niter);n++) 
    {
      mwdebug("iteration : %i / %i \n", n+1,(*Niter));     
      mwdebug("Running InfSup...\n");
      ResolutionInfSup(image,image_work,G,t1,t2,fmovie,Nmask);

      if (average)
	{
	  mwdebug("Running SupInf...\n");
	  ResolutionInfSup(image_save,image_work,G,t2,t1,fmovie,Nmask);	  
	  size = image->ncol*image->nrow;
	  mwdebug("Averaging...\n");
	  for (ptr1=image->gray, ptr2=image_save->gray, j=0; 
	       j<size; ptr1++,ptr2++,j++)
	    *ptr1 = nint(((double) *ptr1 + *ptr2)/2.0);
	  mw_copy_cimage(image,image_save);
	}
    }

  mw_copy_cimage(image,output); 

  /* -------------- Free buffers ------------ */
  
  cfree(G);
  if (image_save) mw_delete_cimage(image_save);
  mw_delete_cimage(image_work);

  return(output);
}