rm libmoduleDfBF

libmoduleDfBF/diffusion_X1.c

-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include "nrutil.h"
-
-#define ISSETBITFLAG(x,b) ((x) & (1 << (b)))
-#define ADD_DIFFUSION    1
-#define ADD_BF_FILTER    2
-
-float linearInterp(float xp, float x0, float y0, float x1, float y1);
-
-void addEffects(int ngx_ima, int ngy_ima, float *arr_ima, float *arr_imc, int bit_flag)
-{
-  int nx, ny, i,j,k,ks;
-  int it,jt,itt,jtt;
-  int diffuidx[26][2],diffuN,ilow,ih,im,dim[3];
-  float diffua[5][5],cdiffu[26],**bfa;
-  double mvar,mcov,tmp,ma,mb,mc;
-  char fname[100];
- 
-  nx = ngx_ima;  //input-image size
-  ny = ngy_ima; 
-
-  //0. init. original image with an input array (arr_ima)
-  //1. Adding diffusion effect.
-  if(ISSETBITFLAG(bit_flag, ADD_DIFFUSION)) 
-  {
-    printf("adding diffusion.....\n");
-    printf("ERR: no diffusion filter ...");
-    exit(0);
-  }
-
-
-  //2. Adding BF effect
-  if(ISSETBITFLAG(bit_flag, ADD_BF_FILTER)) 
-  {
-    printf("Adding BF effect...\n");
-    //setup BF correlation fliter
-    float neX;
-    float neP1 = 50000;
-    float bfaP1[9]={0.9707182, 0.002143905, 0.004131103, 0.001149542, 0.0005501739, 0.0005469659, 0.0003726081, 0.0003795207, 0.0001633302};
-    float neP2 = 10000;
-    float bfaP2[9]={0.9945288, 0.0003041936, 0.0007539311, 0.0002424675, 0.0001226098, 0.00009308617, 0.00008027447, 0.00006309676, 0.00006400052};
-
-    bfa=matrix(-2,2,-2,2);
-
-    // smooth with the BF filter
-    for(i=0;i<nx;i++)for(j=0;j<ny;j++) arr_imc[j+i*ny]=0;
-
-    for(i=0;i<nx;i++)
-    {
-      for(j=0;j<ny;j++)
-      {
-        //rescale BF filter with the local pix value
-        neX = arr_ima[j+i*ny];
-        if(neX >= 10000)
-        {
-          bfa[0][0]=0; //linearInterp(neX, neP1, bfaP1[0], neP2, bfaP2[0]); //0;
-          bfa[0][1]=bfa[0][-1]=linearInterp(neX, neP1, bfaP1[1], neP2, bfaP2[1]); //0.01575;
-          bfa[-1][0]=bfa[1][0]=linearInterp(neX, neP1, bfaP1[2], neP2, bfaP2[2]); //0.00652;
-          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=linearInterp(neX, neP1, bfaP1[3], neP2, bfaP2[3]); //0.00335;
-          bfa[0][-2]=bfa[0][2]=linearInterp(neX, neP1, bfaP1[4], neP2, bfaP2[4]);
-          bfa[-2][0]=bfa[2][0]=linearInterp(neX, neP1, bfaP1[5], neP2, bfaP2[5]); //0.00118;
-          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=linearInterp(neX, neP1, bfaP1[6], neP2, bfaP2[6]);
-          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=linearInterp(neX, neP1, bfaP1[7], neP2, bfaP2[7]); //0.00083;
-          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=linearInterp(neX, neP1, bfaP1[8], neP2, bfaP2[8]); //0.00043;
-        }
-        else
-        {
-          neX=10000;
-          bfa[0][0]=0;
-          bfa[0][1]=bfa[0][-1]=bfaP2[1];
-          bfa[-1][0]=bfa[1][0]=bfaP2[2];
-          bfa[-1][-1]=bfa[1][1]=bfa[-1][1]=bfa[1][-1]=bfaP2[3];
-          bfa[0][-2]=bfa[0][2]=bfaP2[4];
-          bfa[-2][0]=bfa[2][0]=bfaP2[5];
-          bfa[-2][-1]=bfa[-2][1]=bfa[2][1]=bfa[2][-1]=bfaP2[6];
-          bfa[-1][-2]=bfa[1][2]=bfa[-1][2]=bfa[1][-2]=bfaP2[7];
-          bfa[-2][-2]=bfa[-2][2]=bfa[2][-2]=bfa[2][2]=bfaP2[8];
-        }
-        tmp = 0;
-        for(it=-2;it<=2;it++)
-          for(jt=-2;jt<=2;jt++)
-          {
-            bfa[it][jt] = bfa[it][jt]/neX*arr_ima[j+i*ny];
-            tmp += bfa[it][jt];
-          }
-        bfa[0][0]=1.-tmp;
-
-        // assign electrons according to the BF filter bfat
-        for(it=-2;it<=2;it++)
-        {
-          for(jt=-2;jt<=2;jt++)
-          {
-            itt=i+it;
-            jtt=j+jt;
-            if(itt>=0 && jtt>=0 && itt<nx && jtt<ny)
-              //c0[itt][jtt]+=bfa[it][jt]*b[i][j];
-              arr_imc[jtt+itt*ny] += bfa[it][jt]*arr_ima[j+i*ny];
-          }
-        }
-      }
-    }
-    free_matrix(bfa,-2,2,-2,2);
-  }
-  else
-  {
-    for(i=0;i<nx;i++) for(j=0;j<ny;j++) arr_imc[j+i*ny]=arr_ima[j+i*ny]; ////for ADD_BF False
-  }
-}
-
-
-float linearInterp(float xp, float x0, float y0, float x1, float y1)
-{
-  float yp;
-  yp = y0 + ((y1-y0)/(x1-x0)) * (xp - x0);
-  return yp;
-}
-