Codes of Richardson Lucy Algorithm

Reading The Pixel Values:

> library(pixmap)
> lena <- read.pnm("lenna_blur.pgm")
> write.table(lena@grey,"mylenna", quote=FALSE, row.names = FALSE, col.names= FALSE)

C Code:

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

#define H 7
#define PIX 512
#define NUM 512 * 512
#define SIGMA 10

#ifndef max
  #define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

#ifndef min
  #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

void dataread(float *y);
int distance(int i, int j, int k, int l);
void lambda(float *y, float *lold, float *lnew, float *p);
void renew(float *lold,float *lnew);
void spread(float *p);
void datasave(float *l);

int main(void){

  int i;
  float *y, *lold, *lnew, *p, ttt;
  y = malloc(NUM * sizeof(float));
  lold = malloc(NUM * sizeof(float));
  lnew = malloc(NUM * sizeof(float));
  p = malloc((2 * (H + 1) * (H + 1)) * sizeof(float));

  spread(p);
  dataread(y);
  dataread(lold);

  for(i = 0; i < 100; i++){
    lambda(y, lold, lnew, p);
    renew(lold, lnew);
    printf("%d  \n",i);
  }

  datasave(lnew);
  
  free(p);
  free(y);
  free(lnew);
  free(lold);

  return 0;
}


/*READING THE DATA*/

void dataread(float *y){
  int i;
  FILE *fp;
  fp = fopen("mylenna","r");
  for(i = 0; i < NUM; i++){
    fscanf(fp,"%f ",&y[i]);
  }
  return;
}

/*DEFINING A NORM*/

int distance(int i, int j, int k, int l){
  int d;
  d =  (i - k) * (i - k) + (j - l) * (j - l);
  return d;
}




/*LAMBDA VECTOR*/

void lambda(float *y, float *lold, float *lnew, float *p){
  int i1, i2, j1, j2, k1, k2;
  float temp, tmp;

    for(i1 = 0; i1 < PIX; i1++){
      for(i2 = 0; i2 < PIX; i2++){

	tmp = 0;
	for(j1 = max(0, i1 - H + 1); j1 < min(i1 + H - 1, PIX); j1++){
	  for(j2 = max(0, i2 - H + 1); j2 < min(i2 + H - 1, PIX); j2++){

	    temp = 0;

	    for(k1 = max(0, i1 - H + 1); k1 < min(i1 + H - 1, PIX); k1++){
	      for(k2 = max(0, i2 - H + 1); k2 < min(i2 + H - 1, PIX); k2++){
		
		temp = temp + lold[k1 * PIX + k2] * p[distance(j1,j2,k1,k2)];
	      }
	    }
	    
	    tmp = tmp + y[j1 * PIX + j2] * (p[distance(j1,j2,i1,i2)])/temp;
	  }
	}
	lnew[i1 * PIX + i2] = lold[i1 * PIX + i2] * tmp;
      }
    }
  return;
}


/*Point Spread Function*/

void spread(float *p){
  int i, j;
  float temp = 0;
  for(i = -H; i <= H; i++){
    for(j = -H; j <= H; j++){
	p[i * i + j * j] =  exp(-((float)(i * i + j * j))/SIGMA);
	temp += p[i * i + j * j];
    }
  }
  
    
  for(i = -H; i <= H; i++)
    for(j = -H; j <= H; j++)
      p[i * i + j * j] = p[i * i + j * j]/temp;
  return;
}

/*RENEW THE VALUE OF LAMBDA VECTOR*/

void renew(float *lold,float *lnew){
  int i;
  for(i = 0; i < NUM; i++)
    lold[i] = lnew[i];
  return;
}


/*SAVING THE DATA*/

void datasave(float *l){
  int i, j;
  FILE *fp;
  fp = fopen("mylenn","w");

  for(i = 0; i < PIX; i++){
    for(j = 0; j < PIX; j++){
      fprintf(fp,"%f ",l[i * PIX + j]);
    }
    fprintf(fp,"\n");
  }

  return;
}

Map of New Pixels

> plot(pixmapGrey(as.matrix(read.table("mylenn"))))

Interface Between R and C:

#include <R.h>
#include <Rdefines.h>
#include <Rinternals.h>
#include <stdlib.h>

#define H 7
#define PIX 512
#define NUM 512 * 512
#define SIGMA 10

#ifndef max
  #define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

#ifndef min
  #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

int distance(int i, int j, int k, int l);
void lambda(double *y, double *lold, double *lnew, double *p);
void renew(double *lold,double *lnew);
void spread(double *p);


SEXP lenna(SEXP data)
{
  int i;
  double *y, *lold, *lnew, *p;
  SEXP ans;

  ans = PROTECT(NEW_NUMERIC(NUM));
  double *z = NUMERIC_POINTER(ans);

  y = malloc(NUM * sizeof(double));
  lold = malloc(NUM * sizeof(double));
  lnew = malloc(NUM * sizeof(double));
  p = malloc((2 * (H + 1) * (H + 1)) * sizeof(double));

  for(i = 0; i < NUM; i++){
    y[i] = NUMERIC_POINTER(data)[i];
    lold[i] = y[i];
  }

  spread(p);

  for(i = 0; i < 100; i++){
    lambda(y, lold, lnew, p);
    renew(lold, lnew);
    /* printf("%d  \n",i); */
  }

  /* datasave(lnew); */
  for(i = 0; i < NUM; i++)
    z[i] = lnew[i];
  
  free(p);
  free(y);
  free(lnew);
  free(lold);

  UNPROTECT(1);
  return ans;
}


/*DEFINING A NORM*/

int distance(int i, int j, int k, int l){
  int d;
  d =  (i - k) * (i - k) + (j - l) * (j - l);
  return d;
}




/*LAMBDA VECTOR*/

void lambda(double *y, double *lold, double *lnew, double *p){
  int i1, i2, j1, j2, k1, k2;
  double temp, tmp;

    for(i1 = 0; i1 < PIX; i1++){
      for(i2 = 0; i2 < PIX; i2++){

	tmp = 0;
	for(j1 = max(0, i1 - H + 1); j1 < min(i1 + H - 1, PIX); j1++){
	  for(j2 = max(0, i2 - H + 1); j2 < min(i2 + H - 1, PIX); j2++){

	    temp = 0;

	    for(k1 = max(0, i1 - H + 1); k1 < min(i1 + H - 1, PIX); k1++){
	      for(k2 = max(0, i2 - H + 1); k2 < min(i2 + H - 1, PIX); k2++){
		
		temp = temp + lold[k1 * PIX + k2] * p[distance(j1,j2,k1,k2)];
	      }
	    }
	    
	    tmp = tmp + y[j1 * PIX + j2] * (p[distance(j1,j2,i1,i2)])/temp;
	  }
	}
	lnew[i1 * PIX + i2] = lold[i1 * PIX + i2] * tmp;
      }
    }
  return;
}


void spread(double *p){
  int i, j;
  double temp = 0;
  for(i = -H; i <= H; i++){
    for(j = -H; j <= H; j++){
	p[i * i + j * j] =  exp(-((double)(i * i + j * j))/SIGMA);
	temp += p[i * i + j * j];
    }
  }
  
    
  for(i = -H; i <= H; i++)
    for(j = -H; j <= H; j++)
      p[i * i + j * j] = p[i * i + j * j]/temp;
  return;
}

/*RENEW THE VALUE OF LAMBDA VECTOR*/

void renew(double *lold,double *lnew){
  int i;
  for(i = 0; i < NUM; i++)
    lold[i] = lnew[i];
  return;
}