#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;
}