fsepconvol.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- MegaWave2 module  -----------------------------*/
/* mwcommand
 name = {fsepconvol};
 author = {"Lionel Moisan"};
 version = {"1.2"};
 function = {"Convolution with a separable 2D-kernel"};
 usage = {
   'c':width->width           
    "unit mass constant kernel of specified width",
   'g':std->std[0.0,1.0e4]    
    "Gaussian kernel, standart deviation std (precision 1e-3)",
   'b':[b=1]->b  "boundary condition: 0=zero, 1=symmetry (default), 2=torus",
   in->in        "input Fimage",
   out<-out      "output (convolved Fimage)",
 { xker->xker    "use a Fsignal (and its shift) as kernel along x axis",
   yker->yker    "use a Fsignal (and its shift) as kernel along y axis" }
};
*/
/*----------------------------------------------------------------------
 v1.2: new sgauss(), -b option, acceleration, return output (L.Moisan)
----------------------------------------------------------------------*/

/* NOTE: calling this module with in=out is possible */

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

extern void sconst();
extern Fsignal sgauss();


Fimage fsepconvol(in,out,xker,yker,width,std,b)
Fimage in,out;
Fsignal xker,yker;
int *width,*b;
float *std;
{
  Fimage tmp;
  int nx,nx2,ny,ny2,x,y,org,i,iter,s;
  float value,sum,precision;
  int delete_xker;

  /* check options */
  if ((xker?1:0) + (width?1:0) + (std?1:0) != 1)
    mwerror(USAGE,1,"Please specify exactly one kernel");

  if (!xker) {
    xker = mw_new_fsignal();
    delete_xker = 1;
  } else delete_xker = 0;

  nx = in->ncol; nx2 = 2*nx;
  ny = in->nrow; ny2 = 2*ny;
  tmp = mw_change_fimage(NULL,ny,nx);
  out = mw_change_fimage(out,ny,nx);
  
  /*** COMPUTE KERNEL IF NEEDED ***/

  if (width) {
    /* constant kernel */
    value = 1.0/(float)(*width);
    sconst(width,&value,xker);
    xker->shift = -0.5*(float)(*width);
    yker = xker;
  }
  if (std) {
    /* Gaussian kernel */
    precision = 3.;
    sgauss(*std,xker,NULL,&precision);
    yker = xker;
  }

  /* convolution along x axis */
  org = -(int)xker->shift;
  for (y=ny;y--;) for (x=nx;x--;) {
    for (sum=0.0,i=xker->size;i--;) {
      s = x-i+org;
      switch(*b) {
      case 0: 
	if (s>=0 && s<nx) sum += xker->values[i]*in->gray[y*nx+s];
	break;
      case 1: 
	while (s<0) s+=nx2;
	while (s>=nx2) s-=nx2;
	if (s>=nx) s = nx2-1-s;
	sum += xker->values[i]*in->gray[y*nx+s];
	break;
      case 2: 
	while (s<0) s+=nx;
	while (s>=nx) s-=nx;
	sum += xker->values[i]*in->gray[y*nx+s];
	break;
      }
    }
    tmp->gray[y*nx+x] = sum;
  }
  
  /* convolution along y axis */
  org = -(int)yker->shift;
  for (y=ny;y--;) for (x=nx;x--;) {
    for (sum=0.0,i=yker->size;i--;) {
      s = y-i+org;
      switch(*b) {
      case 0: 
	if (s>=0 && s<ny) sum += yker->values[i]*tmp->gray[s*nx+x];
	break;
      case 1: 
	while (s<0) s+=ny2;
	while (s>=ny2) s-=ny2;
	if (s>=ny) s = ny2-1-s;
	sum += yker->values[i]*tmp->gray[s*nx+x];
	break;
      case 2: 
	while (s<0) s+=ny;
	while (s>=ny) s-=ny;
	sum += yker->values[i]*tmp->gray[s*nx+x];
	break;
      }
    }
    out->gray[y*nx+x] = sum;
  }

  /* free memory */
  mw_delete_fimage(tmp);
  if (delete_xker) 
    mw_delete_fsignal(xker);

  return(out);
}