Sourcecode

/**************************************************************************
 ***
 *** PROGRAM: invert_cut.cc
 ***
 *** Program developed for Project 1 CSE 585
 ***
 *** This program can degrade and restore images using various
 *** methods. It was written for the second part of project 1.
 *** 
 *** Author:                Nils Krahnstoever
 *** Date of last revision: 3/3/99
 ***
 **************************************************************************/

/**************************************************************************
 ***
 *** INCLUDES
 ***
 **************************************************************************/

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

// Numerical Recipes Include Files

#define NRANSI
#include "nr.h"
#include "nrutil.h"  
#include "complex.h"

// Personal classes for image and matrix handling

#include <matrix.h>
#include <matrixtools.h> 
#include <image.h>
#include <imagetools.h>

// Personal classes for converting images to matrices and vice versa
#include <converttools.h>

/**************************************************************************
 ***
 *** GLOBAL VARIABLES
 ***
 **************************************************************************/

// Global variables storing the images of differnet stages of the 
// restoration.
Matrix<float> img_re(256,256);  // working matrix
Matrix<float> img_im(256,256);
Matrix<float> ori_re(256,256);  // blurred noise free image
Matrix<float> ori_im(256,256);
Matrix<float> sys_re(256,256);  // system response
Matrix<float> sys_im(256,256);
Matrix<float> res_re(256,256);  // working matrix (often for results)
Matrix<float> res_im(256,256);
Matrix<float> res_abs(256,256); // absolute results
Matrix<float> lap_re(256,256);  // laplacian
Matrix<float> lap_im(256,256);
Matrix<float> noise_fft(256,256);  // fft of noise

// Global variables storing the input arguments.
int mode;        // restoration mode
float value;     // restoration parameter
int flag;        // general flags
int len;         // length of motion blurring or -1 if defocus
float dir;       // direction of motion blurring or defocus sigma
float sigma;     // noise sigma
float from,to;   // repeat restoration varying <value> from <from> to <to>...
int steps;       // ... in <steps> steps

/**************************************************************************
 ***
 *** FUNCTIONS
 ***
 **************************************************************************/

/**************************************************************************
 ***
 *** gaussnoise
 ***
 *** returns a gaussian distributed random variable of mean <mean>
 *** and standard deviation <var>
 ***
 **************************************************************************/
float gaussnoise(float mean, float var)
{
  static int first=1;
  static long seed;

  if(first)
    {
      first=0;
      seed=(long)&seed;
      if(seed>0) seed=-seed;
      ran2(&seed);
    }  

  return(mean+var*gasdev(&seed));
}

/**************************************************************************
 ***
 *** add_noise
 ***
 *** add gaussian distributed noise of mean 0 and stdev of 
 *** <sigma> (global) to the matrix <mat>
 *** 
 **************************************************************************/
void add_noise(Matrix<float> &mat)
{
  static int fl=0;
  int rs=mat.rows();
  int cs=mat.cols();
  int r;
  int c;
  static Matrix<float> re(256,256);
  static Matrix<float> im(256,256);
  static Image *imgptr;

  for(r=0;r<rs;r++)
    {
      for(c=0;c<cs;c++)
        {
          re[r][c]=gaussnoise(0,sigma);
          mat[r][c]+=re[r][c];
        }
    }

  if(fl==0) // the first time this function is called
    {       // create an FFT of the noise
      fl=1;

      im.zero();

      fft(re,im,0);             // fft
      Cabs(noise_fft,im,re);    // abs
      lnp1(noise_fft);          // ln( . + 1 ) => ln(|image|+1)
      
      if(flag & 1) // save image of fft
        {
          imgptr=range_matrix2image(noise_fft);
          imgptr->saveimg("noise_fft.pgm");    
          delete imgptr;
        }
    }
}

/**************************************************************************
 ***
 *** create_defocus
 ***
 *** create a system response function of a blurring process
 *** caused by a gaussian with sigma of <dir> (global)
 *** 
 **************************************************************************/
void create_defocus(Matrix<float> &re,Matrix<float> &im)
{
  int w=256,h=256,x,y,x0,y0;
  float p=acos(-1),sigma,u,v,ua;
  double norm;
  Matrix<float> abs(h,w);
  Matrix<float> re2(h,w),im2(h,w);
  Image *img;

  re.zero();
  im.zero();

  fprintf(stderr,"Creating response (defocus)\n");

  sigma=2*p*p*dir*dir;  // correct sigma in frequency space
  norm=0.0;

  for(x=0;x<w;x++)
    {
          for(y=0;y<h;y++)
            {
              x0=((x+127)%256-127); // origin is at [0][0]
              y0=((y+127)%256-127);
              u=float(x0)/w;
              v=float(y0)/h;
              ua=u*u+v*v;
              re[y][x]= exp(-ua*sigma);
              norm+=re[y][x];
              im[y][x]=0.0;
            }
    }

  // output various images (for debugging)
  if(flag & 1)
    {
      Cabs(abs,re,im);
      lnp1(abs);
      img=range_matrix2image(abs);
      //img=matrix2image(abs);
      img->saveimg("response.pgm");
      delete img;
      re2=re;
      im2=im;
      fft(re2,im2,1);
      Cabs(abs,re2,im2);
      abs.scale(1.0/256/256);
      stat(abs,stderr);
      //img=matrix2image(abs);
      img=range_matrix2image(abs);
      img->saveimg("response_fft.pgm");
      delete img;
    }
}

/**************************************************************************
 ***
 *** create_motion
 ***
 *** create a system response function of a blurring process
 *** caused by motion of length <len> pixel in direction <dir>
 ***
 *** two methods are used:
 ***
 *** if bit 3 of <flag> is set the response function is calculated
 *** analytically
 ***
 *** if bit 3 of <flag> is zero, the convolution mask is created directly
 *** and the response function obtained by FFTing the mask
 *** 
 *** Throughout this project bit 3 of <flag> was set!
 **************************************************************************/
void create_motion(Matrix<float> &re,Matrix<float> &im)
{
  int w=256,h=256,x,y,x0,y0;
  int t;
  float p=acos(-1),dx,dy,u,v,ua,t2,si;
  Matrix<float> abs(h,w);
  Image *img;

  re.zero();
  im.zero();

  fprintf(stderr,"Creating response (motion)\n");

  if(flag & 8)
    {
      dx=cos(dir)*p;
      dy=-sin(dir)*p;
      for(x=0;x<w;x++)
        {
          for(y=0;y<h;y++)
            {
              x0=((x+127)%256-127);
              y0=((y+127)%256-127);
              u=float(x0)/w;
              v=float(y0)/h;
              ua=(dx*u+dy*v)*len;
              if(fabs(ua)<1e-15) 
                {
                  si=1.0;
                }
              else
                {
                  si=sin(ua)/ua;
                }
              re[y][x]=si*cos(-ua);
              im[y][x]=si*sin(-ua);
            }
        }
    }
  else
    {
      for(t=0;t<len;t++)
        {
          t2=t-len/2;
          x0=256+(int)(t2*cos(dir));
          y0=256+(int)(-t2*sin(dir));
          re[y0%256][x0%256]=1.0/len;
        }      
      fft(re,im,0);
    }
  
  // save image (debugging)
  if(flag & 1)
    {
      Cabs(abs,re,im);
      lnp1(abs);
      img=range_matrix2image(abs);
      img->saveimg("response.pgm");
      delete img;
    }
}

/**************************************************************************
 ***
 *** create_response
 ***
 *** create a system response function of a blurring process
 ***
 **************************************************************************/
void create_response(Matrix<float> &re,Matrix<float> &im)
{
  if(len==-1)
    {
      create_defocus(re,im);
    }
  else
    {
      create_motion(re,im);
    }
}

/**************************************************************************
 ***
 *** degrade_image
 ***
 *** degrade <img> with system response function <sys>
 ***
 **************************************************************************/
Matrix<float> *degrade_image(Matrix<float> &img,Matrix<float> &sys_re,Matrix<float> &sys_im)
{
  int r,c;

  Matrix<float> *mat=new Matrix<float>(256,256);
  Matrix<float> img_im(256,256);
  Matrix<float> img_re(256,256);

  fcomplex a,b,s;
  img_re=img;
  img_im.zero();

  fprintf(stderr,"norm(image)=%g   ",sum(img_re));

  // <img> -> FFT domain
  fft(img_re,img_im,0);

  // FFT(<img>)*<sys>
  r=img_re.rows();
  c=img_re.cols();
  for(r=0;r<256;r++)
    {
      for(c=0;c<256;c++)
        {
          a.r=img_re[r][c];
          a.i=img_im[r][c];
          b.r=sys_re[r][c];
          b.i=sys_im[r][c];
          s=Cmul(a,b);
          img_re[r][c]=s.r;
          img_im[r][c]=s.i;
        }
    }

  // FFT(<blurred image>) -> spatial domain
  fft(img_re,img_im,1);

  // make real
  Cabs(*mat,img_re,img_im);

  // fix scaling
  mat->scale(1.0/256/256);

  // did the norm change much? (debugging)
  fprintf(stderr,"norm(deg)=%g\n",sum(*mat));

  // return degraded image
  return(mat);
}

/**************************************************************************
 ***
 *** laplacian
 ***
 *** if mode==4 (constr. least square) we need a laplacian in freq. space
 *** this function creates one
 ***
 **************************************************************************/
void laplacian(Matrix<float> &lap_re,Matrix<float> &lap_im)
{
  Matrix<float> lap_abs(256,256);
  Image *imgptr;

  lap_re.zero();
  lap_im.zero();

  //  0 -1  0
  // -1  4 -1
  //  0 -1  0

  lap_re[0][0]=0;
  lap_re[0][1]=-1;
  lap_re[0][2]=0;
  lap_re[1][0]=-1;
  lap_re[1][1]=4;
  lap_re[1][2]=-1;
  lap_re[2][0]=0;
  lap_re[2][1]=-1;
  lap_re[2][2]=0;

  // FFT
  fft(lap_re,lap_im,0);

  // save (debugging)
  if(flag & 1)
    {
      Cabs(lap_abs,lap_re,lap_im);
      lnp1(lap_abs);
      imgptr=range_matrix2image(lap_abs);
      imgptr->saveimg("laplacian.pgm");    
      delete imgptr;
    }
}

/**************************************************************************
 ***
 *** invert
 ***
 *** invert image <img> with respect to system response <sys>
 ***
 *** <mode> (global variable) denotes the method used
 ***
 **************************************************************************/
Matrix<float> *invert(Matrix<float> &img_re,Matrix<float> &img_im,
                      Matrix<float> &sys_re,Matrix<float> &sys_im)
{
  fcomplex a,b,s;
  float dx,d,p,v,n,f;
  int r,c;
  Matrix<float> res_re(256,256);  // result
  Matrix<float> res_im(256,256);
  Matrix<float> *res;             // result to be returned

  res=new Matrix<float>(256,256); // create matrix

  if(mode==2) v=value*value;      // mode2 cutoff frequency^2

  for(r=0;r<256;r++)              // for each pixel
    {
      for(c=0;c<256;c++)
        {
          a.r=img_re[r][c];       // a -> G(u,v)
          a.i=img_im[r][c];
          b.r=sys_re[r][c];       // b -> H(u,v)
          b.i=sys_im[r][c];
          // G/H and the variouse flavor of this are stored in s
          switch(mode)
            {
            case 1:               // invert ignoring values less than <value>
              {
                if(Cabs(b)<value)
                  {
                    s.r=a.r;
                    s.i=a.i;
                  }
                else
                  s=Cdiv(a,b);
                break;
              }
            case 2:               // invert with frequency cut-off
              {
                dx=((r+127)%256-127);
                dx*=dx;
                d=((c+127)%256-127);
                d*=d;
                d+=dx;
                if(d>v || Cabs(b)<0.05)  // still ignore small values
                  {
                    s.r=a.r;
                    s.i=a.i;
                  }
                else
                  s=Cdiv(a,b);
                break;
              }
            case 3:              // add constant (Wiener with constant)
              {
                d=Cabs(b);
                d=1/(d*d+value);
                b.r*=d;
                b.i*=-d;
                s=Cmul(a,b);
                break;
              }
            case 4:              // constr. least square using
                                 // laplacian <lap> 
              {
                d=Cabs(b);
                s.r=lap_re[r][c];
                s.i=lap_im[r][c];
                p=Cabs(s);
                d=1/(d*d+value*p*p);
                b.r*=d;
                b.i*=-d;
                s=Cmul(a,b);
                break;
              }     
            case 6:              // power Spectrum Filter
              {
                d=Cabs(b);
                n=exp(noise_fft[r][c])-1;  // noise est.
                // original image (quite unrealistic)
                f=ori_re[r][c]*ori_re[r][c]+ori_im[r][c]*ori_im[r][c];
                d=sqrt(f/(d*d*f+value*n*n));
                s.r=a.r*d;
                s.i=a.i*d;
                break;
              }     
            case 7:              // parametric Wiener Filter
              {
                d=Cabs(b);
                n=exp(noise_fft[r][c])-1;  // noise est.
                // original image (quite unrealistic)
                f=ori_re[r][c]*ori_re[r][c]+ori_im[r][c]*ori_im[r][c];
                d=1/(d*d+value*n*n/f);
                b.r=b.r*d;
                b.i=-b.i*d;
                s=Cmul(a,b);
                break;
              }     
            default:             // just divide
              s=Cdiv(a,b);
            }
          res_re[r][c]=s.r;
          res_im[r][c]=s.i;
        }
    }

  // FFT back
  fft(res_re,res_im,1);

  // make real
  Cabs(*res,res_re,res_im);

  // rescale
  res->scale(1.0/256/256);

  // return result
  return res;
}

/**************************************************************************
 ***
 *** rms_diff
 ***
 *** calculates the root mean square of difference between the original
 *** image and the <mat>
 ***
 **************************************************************************/
float rms_diff(Image &img,Matrix<float> &mat)
{
  float rms=0.0,d;
  int i,x,y;

  i=0;
  for(y=0;y<256;y++)
    { 
      for(x=0;x<256;x++)
        {
          d=mat[y][x]-(float)img.data(i);
          i++;
          d=d*d;
          rms+=d;
        }
    }
  rms=sqrt(rms/256/256);
  return rms;
}

/**************************************************************************
 ***
 *** var_diff
 ***
 *** calculates the variance of the differences
 *** this measure makes it possible to assure that a high rms error is
 *** not caused by a systematic offset of pixel intensity values in the result
 ***
 **************************************************************************/
float var_diff(Image &img,Matrix<float> &mat)
{
  float mean=0.0,var=0.0,d;
  int i,x,y;
  
  i=0;
  for(y=0;y<256;y++)
    { 
      for(x=0;x<256;x++)
        {
          d=(mat[y][x])-(float)img.data(i);
          i++;
          mean+=d;
        }
    }

  mean=mean/256/256;

  i=0;
  for(y=0;y<256;y++)
    { 
      for(x=0;x<256;x++)
        {
          d=mean-((mat[y][x])-(float)img.data(i));
          i++;
          d*=d;
          var+=d;
        }
    }
  var=var/256/256;
  return sqrt(var);
}


/**************************************************************************
 ***
 *** MAIN 
 ***
 **************************************************************************/
void main(int argc,char *argv[])
{
  float nsq;

  Image img(argv[1]);
  Image *imgptr;

  Matrix<float> *res=NULL;
  Matrix<float> *mat;
  Matrix<float> orig(256,256);

  float phi,diff;
  int step=0,noise,nopt;
  char name[80];
  int vary_noise=0;

  // 11 arguments !!!
  if (argc!=11)
    {
      fprintf(stderr,"usage: invert_cut <image [PGM]> <result PGM> <mode> <value> <length> <dir/sigma> <flag> <noise ampl.> <to> <nval>\n");
      fprintf(stderr,"flag : 1 - save respone, laplacian etc...\n");
      fprintf(stderr,"       2 - optimize least square gamma (mode 4)\n");
      fprintf(stderr,"       4 - don't create motion image <image [PGM]> is already blurred\n");
      fprintf(stderr,"       8 - calculate motion response directly (rather than using FFT)\n");
      fprintf(stderr,"       16 - save animation\n");
      fprintf(stderr,"modes: 0 - don't treat zeros\n");
      fprintf(stderr,"       1 - ignore values less than <value>\n");
      fprintf(stderr,"       2 - ignore freq greater <value>\n");
      fprintf(stderr,"       3 - Wiener filtering with constant <value>\n");
      fprintf(stderr,"       4 - constr. least squared with gamma of <value>\n");
      fprintf(stderr,"       5 - only add noise\n");
      fprintf(stderr,"       6 - power spectrum filter\n");
      fprintf(stderr,"       7 - Wiener Filter\n");
      fprintf(stderr,"<len>          = -1 => add defocus image\n");
      fprintf(stderr,"<noise ampl.> <=  0 => no noise\n");
      exit(1);
    }
  printf("----\n");

  // get input variables

  mode=atoi(argv[3]);
  value=atof(argv[4]);
  from=value;
  len=atoi(argv[5]);
  dir=atof(argv[6]);
  to=atof(argv[9]);
  steps=atoi(argv[10]);
  if(len!=-1) dir*=2.0*acos(-1)/360.0;   // deg to rad
  if(argc>=9) sigma=atof(argv[8]);
  flag=0;
  if(argc>=8) flag=atoi(argv[7]);

  // create response function of blurring process

  create_response(sys_re,sys_im);

  // convert input image to matrix format

  mat=image2matrix(img);
 
  // print statistics about image

  stat(*mat,stderr);

  // obtain fft of image and store in ori_re, ori_im
  ori_re=*mat;
  ori_im.zero();
  fft(ori_re,ori_im,0);

  // save resulting fft if desired
  if(flag & 1)
    {
      Cabs(res_abs,ori_re,ori_im);
      lnp1(res_abs);
      imgptr=range_matrix2image(res_abs);
      imgptr->saveimg("orig_fft.pgm");    
      delete imgptr;
    }

  // optain blurred version of image <img> if bit 2 not set
  // by multiplying with the response function
  img_re=*mat;
  if(!(flag & 4))
    {
      fprintf(stderr,"Degrading image!\n");
      if(mode!=5)
        {
          delete mat;
          mat=degrade_image(img_re,sys_re,sys_im); // degrade
        }
    }
  else
    {
      fprintf(stderr,"Not degrading image!\n");
    }

  if(flag & 1)
    {
      Matrix<float> m(256,256);
      m=*mat;
      imgptr=matrix2image(m);
      imgptr->saveimg("motion.pgm");     // save degraded image
      if(sigma>0) add_noise(m);          // add noise
      delete(imgptr);
      imgptr=matrix2image(m);
      imgptr->saveimg("noise.pgm");      // save degraded noisy image
      delete(imgptr);
    }      

  orig=*mat;                        // keep blurred (noise free) image

  // calc fft of blurred image

  img_re=*mat;
  img_im.zero();
  fft(img_re,img_im,0);
  
  if(flag & 1)
    {
      Cabs(res_abs,img_re,img_im);
      lnp1(res_abs);
      imgptr=range_matrix2image(res_abs);
      imgptr->saveimg("image_fft.pgm");  // save fft of blurred image
      delete imgptr;
    }

  // calculate laplacian operator
  if(mode==4)
    {
      laplacian(lap_re,lap_im);
    }
  
  if(sigma<0) vary_noise=1;     // repeat for varying noise levels

  for(noise=0;noise<6;noise++)  // different noise levels (not used in proj.)
    {

      // NOISE CONTROL
      if(vary_noise==1)
        {
          if(noise==0) 
            sigma=0;
          else
            sigma=pow(2.0,(float)noise);  // set noise level
        }
      else
        {
          noise=7;              // don't loop anymore
        }

      // VARY PARAMETERS OF INVERSION
      for(step=0;step<steps;step++)  // change parameter of invert
        {
          // SET <VALUE>
          if(steps<2) 
            value=from;
          else
            value=from+(to-from)*step/(steps-1);

          // OPTIMIZE OF DESIRED
          nopt=0;diff=1e10;
          while(abs(diff)>10 && nopt<10)
                {
                  nopt++;
                  img_re=orig;        // the original blurred image
                  img_im.zero();
                  
                  nsq=(256-1)*(256-1)*sigma*sigma;
                  if(sigma>0) add_noise(img_re);   // add noise

                  fft(img_re,img_im,0);            // get FFT

                  if(res!=NULL) delete res;        // delete old result

                  res=invert(img_re,img_im,sys_re,sys_im);  // invert

                  if(mode==4 && (flag & 2))        // constr. least square?
                    {
                      delete mat;
                      mat=degrade_image(*res,sys_re,sys_im);  // degrade
                      phi=normofdiff(*mat,orig);              // compare
                      diff=(phi-nsq)/256/256;                 // calc diff
                      fprintf(stderr,"OPT %d: diff=%g phi=%g nsq=%g\n",nopt,diff,phi,nsq);
                    }
                  diff=0;               
                }

          // output some quantitative results
          printf("method=%s, par=%2d,  mode=%d,  value=%6.4f,  noise=%02d,  rms error=%8.2f,  variance of diff=%8.2f\n",len < 0 ? "g" : "m",len < 0 ? (int) dir : len,mode,value,(int)sigma,rms_diff(img,*res),var_diff(img,*res));

          // save restored image (parameters part of the name)
          if(flag & 16)
            {
              sprintf(name,"%s_par%03d_mod%d_val%5.4f_noi%02d.pgm"
                      ,len < 0 ? "g" : "m"
                      ,len < 0 ? (int) dir : len
                      ,mode,value,(int)sigma);
              imgptr=mean_matrix2image(*res);
              //imgptr=matrix2image(*res);
              imgptr->saveimg(name);              
              delete imgptr;
            }
        }
    }
  
  // output some statistics
  stat(*res,stderr);

  // transfer to PGM
  imgptr=mean_matrix2image(*res);

  // save final result
  imgptr->saveimg(argv[2]);    

  // delete old images and matrices
  delete imgptr;
  delete res;
}