#include "module.c"
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>


const char *CHANNEL_NAMES[] = {"EVENT_NUM", "EVENT_I", "EVENT_J", "EVENT_CHARGE", "EVENT_ETA_X", "EVENT_ETA_Y", "EVENT_I_INTERP", "EVENT_J_INTERP"};


const int MAX_NUM_EVENTS= 100; // max number of events per frame
const int EVENT_CHANNELS = 7;


struct events_double func_ph_1d_double( double *frame, int i_dim, int j_dim,  double *th_m);
//struct decl.
struct events_double {
  double    **evnt_ijc;
  int evnt_num;
};



//def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None):
PyObject *process(PyObject *self, PyObject *args)
{
    PyArrayObject *image;
    long pulse_id;
    PyObject /*double*/ *timestamp;
    long seconds, nanos;
    PyArrayObject *x_axis;
    PyArrayObject *y_axis;
    PyObject *pars;
    PyObject *bsdata;

    //if (!PyArg_ParseTuple(args, "OldOOO|O", &image, &pulse_id, &timestamp, &x_axis, &y_axis, &pars, &bsdata))
    if ( !PyArg_ParseTuple(args, "OlOOOO|O", &image, &pulse_id, &timestamp, &x_axis, &y_axis, &pars, &bsdata) ||
         !PyArg_ParseTuple(timestamp, "ll", &seconds, &nanos) )
        return NULL;

    if (pulse_id < 0) {
        PyErr_SetString(moduleErr, "Invalid Pulse ID");
        return NULL;
    }

    //Acessing image
    int element_size = image->descr->elsize;
    int dims = image->nd;
    int size_x = image->dimensions[1];
    int size_y = image->dimensions[0];
    unsigned short* img_data = (unsigned short*)image->data;


    int i,j,l;
    int i_dim=size_y;
    int j_dim=size_x;    

    double *threshold = malloc(i_dim*j_dim*sizeof(double));
    for(i=0; i<i_dim;i++) {
      for(j=0;j<j_dim;j++) {
              threshold[i*j_dim+j]=10.0;
          }
       }
    //background (all matrices are indexed in 1d)
    double *background = malloc(i_dim*j_dim*sizeof(double));
    for(i=0; i<i_dim;i++) {
      for(j=0;j<j_dim;j++) {
              background[i*j_dim+j]=210.0;
          }
       }


    double *frameBKsub = malloc(i_dim*j_dim*sizeof(double));
    for(i=0; i<i_dim;i++) {
      for(j=0;j<j_dim;j++) {
              frameBKsub[i*j_dim+j]=(double)img_data[i*j_dim+j] - background[i*j_dim+j];
          }
       }
    struct events_double evns1d = func_ph_1d_double( frameBKsub, i_dim, j_dim, threshold);

    //Create return dictionary
    PyObject *ret = PyDict_New();
    PyObject* camera_name = PyDict_GetItemString(pars, "camera_name");
    const char * camera_name_str = PyUnicode_AsUTF8(camera_name);
    char channel_name[200];        
    //PyDict_SetItemString(ret, "camera_name", camera_name);
    sprintf(channel_name, "%s:%s",camera_name_str,  CHANNEL_NAMES[0]);
    PyDict_SetItemString(ret, channel_name, PyLong_FromLong(evns1d.evnt_num));
    for (int i=0; i<EVENT_CHANNELS; i++){
        sprintf(channel_name, "%s:%s",camera_name_str,  CHANNEL_NAMES[i+1]);
        npy_intp arr_dims[1] = {MAX_NUM_EVENTS};
        PyObject *parr =PyArray_SimpleNewFromData(1, arr_dims, NPY_FLOAT64, evns1d.evnt_ijc);
        PyDict_SetItemString(ret, channel_name, parr);
     }
    for (l = 0; l < EVENT_CHANNELS; l++) {
       free(evns1d.evnt_ijc[l]);
     }
     free(evns1d.evnt_ijc);
     free(frameBKsub);
     free(background);
     free(threshold);
    return ret;
}

 struct events_double func_ph_1d_double( double *frame, int i_dim, int j_dim,  double *th_m)
 {
   //int th= 50;
   int i, j, l=0, m=0, n=0, evt_m_i=0, evt_i=0 ;
   int dist_i, dist_j;
   double charge_evt;
   int th;

   double **evt_p= (double**)malloc(EVENT_CHANNELS * sizeof(double*));
   for (l = 0; l < EVENT_CHANNELS; l++) {
       evt_p[l] = (double*)malloc(MAX_NUM_EVENTS * sizeof(double));}

   /*Counter variables for the loop*/
   double charge = 0;
   double eta_x =  0;
   double eta_y =  0;
   double i_interp = 0;
   double j_interp = 0;
   double C[]= { 4.37805097e-03, 2.43266401e-01, -7.81479328e+00, 7.70533057e+01, -2.37906124e+02, 3.42988113e+02, -2.35979731e+02, 6.24074562e+01};
   double D[]= {-2.05180773e-04, 4.35314696e-01, -1.32557223e+01, 1.19990373e+02, -3.80696758e+02, 5.81427096e+02, -4.30789104e+02, 1.23894761e+02};

   for(i=0; i<i_dim-1; i++) {      
          if (evt_i>=MAX_NUM_EVENTS){
               break;
          }
      for(j=0;j<j_dim-1;j++) {
          if (evt_i>=MAX_NUM_EVENTS){
               break;
          }


          // 2x2 version
          charge = frame[i*j_dim+j]+frame[(i+1)*j_dim+j] + frame[i*j_dim+(j+1)]+frame[(i+1)*j_dim+j+1];

          //pixel by pixel threshold
          th = th_m[i*j_dim +j];

          //check if charge above threshold
          if(charge>th) {
                 eta_x    = (frame[(i+1)*j_dim + j    ]+frame[(i+1)*j_dim +  (j+1)])/charge;
                 eta_y    = (frame[    i*j_dim + (j+1)]+frame[(i+1)*j_dim +  (j+1)])/charge;
                 i_interp =  i + (C[0] + C[1]*eta_x+ C[2]*pow(eta_x,2) + C[3]*pow(eta_x,3) + C[4]*pow(eta_x,4)+ C[5]*pow(eta_x,5) + C[6]*pow(eta_x,6) + C[7]*pow(eta_x,7));
                 j_interp =  j + (D[0] + D[1]*eta_y+ D[2]*pow(eta_y,2) + D[3]*pow(eta_y,3) + D[4]*pow(eta_y,4)+ D[5]*pow(eta_y,5) + D[6]*pow(eta_y,6) + D[7]*pow(eta_y,7));

                 // 1st case: first event
                 if(evt_i==0){
                 evt_p[0][evt_i] = i;
                 evt_p[1][evt_i] = j;
                 evt_p[2][evt_i] = charge;
                 evt_p[3][evt_i] = eta_x;
                 evt_p[4][evt_i] = eta_y;
                 evt_p[5][evt_i] = i_interp; //
                 evt_p[6][evt_i] = j_interp; //
                 evt_i++;

                 } else {
                 // 2nd case: not 1st event. we check if it is a neighbourg of the previos events and if charge is larger
                 n=0;
                 evt_m_i = evt_i;
                 for(m=0; m<evt_m_i; m++) {
                        dist_i     = abs(evt_p[0][evt_i-1-m] - i); //fix here
                        dist_j     = abs(evt_p[1][evt_i-1-m] - j);
                        charge_evt =     evt_p[2][evt_i-1-m];

                        if( dist_i< 2 && dist_j<2) {
                            if(charge_evt < charge) {
                             //2nd case: not 1st event, but neigboor of previos event, this event has more charge
                             evt_p[0][evt_i-1-m] = i;
                             evt_p[1][evt_i-1-m] = j;
                             evt_p[2][evt_i-1-m] = charge;
                             evt_p[3][evt_i-1-m] = eta_x;
                             evt_p[4][evt_i-1-m] = eta_y;
                             evt_p[5][evt_i-1-m] = i_interp;
                             evt_p[6][evt_i-1-m] = j_interp;
                            } else {
                                //3d case not 1st event, but neigboor of previos event, this event has less charge
                            }

                        } else {
                             // not a neighbor of the m event

                             n++;

                             //now we check if it is not a neighborg of any previous event
                             if((n)==evt_i) {

                             evt_p[0][evt_i] = i;
                             evt_p[1][evt_i] = j;
                             evt_p[2][evt_i] = charge;
                             evt_p[3][evt_i] = eta_x;
                             evt_p[4][evt_i] = eta_y;
                             evt_p[5][evt_i] = i_interp; //
                             evt_p[6][evt_i] = j_interp; //
                             evt_i++;
                             }
                        }


                   }
            }
        }
      }
    }
    struct events_double events_i_d;
    events_i_d.evnt_ijc = evt_p;
    events_i_d.evnt_num = evt_i;
    return events_i_d;
 }