/**************************************************************/
/*                                                            */
/*                        tectum_library.C                    */
/*                                                            */
/**************************************************************/

# include "nsl_include.h"

/* This file involves all the functions concerning the ganglion
   cells; R2, R3 and R4. It gets the retinal input, the black box,
   according to the stimuli type and size. */

/* This is the matrix storing the 'black box' values for the 
   ganglion cells according to the type of stimulus present. */

/*                                  size                        type  */
/*                         2     4     8     16    32                 */
static
float dat[3][3][5]= { { { 25.0, 24.0, 29.0, 28.0, 27.0 },     /* R2 a */
                        { 25.0, 35.0, 25.0, 10.0,  4.0 },     /*    b */
                        { 25.0, 36.0, 23.0,  9.0,  3.0 } },   /*    c */
                      { { 16.0, 20.0, 19.0, 18.0, 20.0 },     /* R3 a */
                        { 16.0, 25.0, 35.0, 17.0, 10.0 },     /*    b */
                        { 16.0, 25.0, 34.0, 20.0,  9.0 } },   /*    c */
                      { { 10.0, 11.0, 10.0, 12.0, 11.0 },     /* R4 a */
                        { 10.0, 16.0, 17.0, 17.0, 15.0 },     /*    b */
                        { 10.0, 17.0, 20.0, 24.0, 17.0 } } }; /*    c */

/* 'get_size' calculates the size to be passed as argument to the
   'black box'. */

int get_size(int d)
{
	int size;

	switch (d)
	{
		case 2  : size = 1;
			  break;
		case 4  : size = 2;
			  break;
		case 8  : size = 3;
			  break;
		case 16 : size = 4;
			  break;
		case 32 : size = 5;
			  break;
		default : size = 0;
			  break;
	}
	return(size);
}

/* The 'black_box' contains the values to be used by the ganglion cells. 

	val0 = vk(V**d) => val = (vk/N)(v**d)

  'd' is the power constant,
  'N' is the normalization factor,
  'vk' is a constant,
  'val0' is the black box frequency response,
  'V' is the speed used to obtain the black box values,
  'v is the stimulus new speed,
  'val' is the resulting frequency response.
*/

float black_box(int n,float dx,float dy,float vx,float d)
{
	static float N=25.0, V=7.6;
	int type,size,idx,idy;
	int i,j,k;
	float val,val0,v,vk;

	idx = nsl_rint(dx);
	idy = nsl_rint(dy);
	v = vx;	/* only x direction (vx) */
	type = 0;
	if (idx == 2 && idy != 2) 	/* WORM */
	{
		type = 1;
		size = get_size(idy);
	}
	if (idx != 2 && idy == 2) 	/* ANTIWORM */
	{
		type = 2;
		size = get_size(idx);
	}
	if (idx == idy)			/* BOX */
	{
		type = 3;
		size = get_size(idx);
	}
	i = n-2;
	j = type-1;
	k = size-1;
	if (i>=0 && i<=2 && j>=0 && j<=2 && k>=0 && k<=4)
	{
		val0 = dat[i][j][k];
        vk = val0/pow(V,d);	/* pow(V,d) => V**d, <math.h> */ 	
        val = vk*pow(v,d)/N;		
	}
    return(val);
}

// 'ganglion_cell_def' is computed at the beginning to get the right 
// 'black box' values for the different ganglion cells.

void ganglion_cell_def(nsl_input_matrix& r,int type,float d,float erf,float val)
{
	float new_val,dx,dy,vx;
	nsl_list* list;
	nsl_stim* obj;

	if ((list = r.get_data_list()) == NULL)
		cmd_error("ganglion_cell_def: no input patterns");
	else
	{
	   list->wind_up(); 
	   while (obj = (nsl_stim*) list->next())
	   {
		dx = ((nsl_block_stim_matrix*) obj)->get_dx();
		dy = ((nsl_block_stim_matrix*) obj)->get_dy();

		vx = ((nsl_stim_matrix*) obj)->get_vx();
//		vy = ((nsl_stim_matrix*) obj)->get_vy();

		if (vx == 0.0)
		   new_val = val;
		else
		   new_val = black_box(type,dx,dy,vx,d);

		((nsl_block_stim_matrix*) obj)->set_val(new_val);
		((nsl_block_stim_matrix*) obj)->set_dx(dx + erf);
		((nsl_block_stim_matrix*) obj)->set_dy(dy + erf);
	   }
	}
}

float dia_chek(float x,float y,float r)
{
        return (x*x+y*y < r*r  ? 1.0 : 0.0);      
}

/* The 'erf_func' function calculates overlapping between the cells
   and the stimulus according to its current position. */

/*

    1		2	    3
A		|		|
  ----------------------
B		|   *	|
  ----------------------
C		|		|

*/

void erf_func(nsl_input_matrix& in,nsl_input_matrix& r,
	int type,float d,float erf,float val)
{                    
	nsl_list* list;

	float tx1,ty1,tx0,ty0,txc,tyc;
	float dx,dy;
    float m,n;
    float ratio;        // val ratio (0.0 to 1.0)
	float vx;
	float new_val;
	nsl_stim* obj;

	int imax = r.get_imax();
	int jmax = r.get_jmax();
	int i0 = r.get_i0();
	int j0 = r.get_j0();
	int i1 = r.get_i1();
	int j1 = r.get_j1();

	float mdx = r.get_dx();
	float mdy = r.get_dy();
	int xz = r.get_xz();
	int yz = r.get_yz();

    float mdx2 = mdx/2; 	// midx interval
    float mdy2 = mdy/2; 	// midy interval
    float rf = erf/2;   	// erf radius
    float rfx = rf/mdx;   	// x-rf radius
    float rfy = rf/mdy;   	// y-rf radius
    float area = mdx*mdy;	// interval area

	static int cnt = 0;
	int fg = 0;

	nsl_input_layer* input_layer = in.get_input_layer();

	if (input_layer != NULL && (list = input_layer->get_stim_list()) != NULL)
	{
	   list->wind_up();
	   while (obj = (nsl_stim*) list->next())
	   {
	    	if (obj->get_able_fg() == 0)
				continue;

			dx = ((nsl_block_stim_matrix*) obj)->get_dx(); // obj width
			dy = ((nsl_block_stim_matrix*) obj)->get_dy(); // obj height

			vx = ((nsl_stim_matrix*) obj)->get_vx();
//			vy = ((nsl_stim_matrix*) obj)->get_vy();

			if (vx == 0.0)
				new_val = val;
			else
				new_val = black_box(type,dx,dy,vx,d);

			txc = ((nsl_stim_matrix*) obj)->get_txc();
			tyc = ((nsl_stim_matrix*) obj)->get_tyc();
			tx0 = ((nsl_stim_matrix*) obj)->get_tx0();
			ty0 = ((nsl_stim_matrix*) obj)->get_ty0();
			tx1 = tx0 + dx;
			ty1 = ty0 + dy;

        	for (int i = i0; i <= i1; i++)
	            for (int j = j0; j <= j1; j++)
	            {   // x,y crossed
	                m = mdx*(j - xz); // current x elem location
	                n = mdy*(i - yz); // current y elem location
	                ratio = 0.0;
	                if ((tx0<m && m<=tx1) && (ty0<n && n<=ty1)) //(tx0<=m && m<=tx1) && (ty0<=n && n<=ty1)) // 
					{
						ratio = 1.0; // 2B
//						if (fg == 1) NSLoutput("\n2 B");
					}
	                else if ((tx0<m && m<=tx1) && // ((tx0<=m && m<=tx1)
	                  ((ty1<n && n-rf<ty1) || (n<ty0 && ty0<n+rf)))
					{
	                    ratio = 1.0; // 2C || 2A
//						if (fg == 1) NSLoutput("\n2 A/C");
					}
	                else if ((ty0<n && n<ty1) && // (ty0<=n && n<=ty1)
	                  ((tx1<m && m-rf<tx1) || (m<tx0 && tx0<m+rf)))
					{
	                    ratio = 1.0; // 3B || 1B
//						if (fg == 1) NSLoutput("\n1/3 B");
					}
	                else if ((tx1<m && ty1<n) || (tx1<m && n<ty0) 
						|| (m<tx0 && n<ty0) || (m<tx0 && ty1<n))
				// 3C || 3A || 1A || 1C
					{
	                    ratio = dia_chek(fabs(m-txc),fabs(n-tyc),rf);
//						if (fg == 1) NSLoutput("\n1/3 A/C");
					}

// fractional value according to overlapping
/*	                if (m-mdx2>=px0 && m+mdx2<=px1)
	                {
	                    if (n-mdy2>=py0 && n+mdy2<=py1) 
	                        r = 1.0;                  // totally covered
	                    else if (py1<n+mdy2 && n-mdy2>=py0) 
	                        r = (py1-(n-mdy2))/mdy;	  // partial upper area
	                    else if (py0<n-mdy2 && n+mdy2<=py1) 
	                        r = ((n+mdy2)-py0)/mdy;   // partial lower area
	                    else if (py1<mdy2+n && py0<mdy2-n)
	                        r = (py1-py0)/mdy;        // thin area
	                }
	                else if (n-mdy2>=py0 && n+mdy2<=py1)
	                {
	                    if (m-mdx2>=px0 && m+mdx2<=px1)
	                        r = 1.0;                  // totally covered
	                    else if (px1<mdx2+m && m-mdx2>=px0)
	                        r = (px1-(m-mdx2))/mdx;   // partial upper area
	                    else if (fabs(m-px0)<mdx2 && m+mdx2<=px1)
	                        r = ((m+mdx2)-px0)/mdx;   // partial lower area
	                    else if (fabs(px1-m)<mdx2 && fabs(m-px0)<mdx2)
	                        r = (px1-px0)/mdx;        // thin area
	                }
	                else if (fabs(px1-m)<mdx2 && fabs(py1-n)<mdy2
	                        && m-mdx2>=px0 && n-mdy2>=py0)
	                    r = (px1-(m-mdx2))*(py1-(n-mdy2));// top right corner
	                else if (fabs(px1-m)<mdx2 && fabs(py1-n)<mdy2
	                        && m-mdx2>=px0 && n-mdy2>=py0)
	                    r = (px1-(m-mdx2))*(py1-(n-mdy2));// top right corner
	                else if (fabs(px1-m)<mdx2 && fabs(py1-n)<mdy2
	                        && m-mdx2>=px0 && n-mdy2>=py0)
	                    r = (px1-(m-mdx2))*(py1-(n-mdy2));// top right corner
	                else if (fabs(px1-m)<mdx2 && fabs(py1-n)<mdy2
	                        && m-mdx2>=px0 && n-mdy2>=py0)
	                    r = (px1-(m-mdx2))*(py1-(n-mdy2));// top right corner
*/
					if (i == 5 && j == 1 && ratio == 1.0)
						;
	                r[i][j] = ratio*new_val;
	            }
	   }
	}
	else
	   cmd_out("NULL 'stim_list'");

	cnt++;
}

/*
int_erf_func(mat,imax,jmax,i0,i1,j0,j1,erf,val)
float **mat;
int imax,jmax,i0,i1,j0,j1;
float erf,val;
{                
        int i,j,r;
        m,n,rf;

        dx = dx/2;
        dy = dy/2;
        rf = erf/2;
        for (i=i0; i<=i1; i++)
           for (j=j0; j<=j1; j++)
           {  // crossed x,y
              m=dxz*(j - xz);
              n=dyz*(i - yz);
              r=0;
              if ((m>=xc-dx && m<=xc+dx) && (n>=yc-dy && n<=yc+dy))
                    r=1;
              if ((m>=xc-dx && m<=xc+dx) && 
                  ((n>yc+dy && n-rf<yc+dy) || (n<yc-dy && n+rf>yc-dy)))
                    r=1;
              if ((n>=yc-dy && n<=yc+dy) &&
                  ((m>xc+dx && m-rf<xc+dx) || (m<xc-dx && m+rf>xc-dx)))
                    r=1;
              if ((m>xc+dx && n>yc+dy) || (m>xc+dx && n<yc-dy) ||
                  (m<xc-dx && n<yc-dy) || (m<xc-dx && n>yc+dy))
                    r=dia_chek(fabs(m-xc)-dx,fabs(n-yc)-dy,rf);
              if (r==1)
                 mat[i][j] = val;
           }
}
*/