#include "stdafx.h"


#using <system.drawing.dll>

using namespace System;
using namespace System::Drawing;
using namespace System::Drawing::Imaging;
//using namespace System::Drawing::Internal;
using namespace std::tr1;


extern  simula_Poisson_noise poi;
extern all_trash AT; 
extern im_GL   Image_0, Image_1,	Image_2; // (2d+t)-input RGBA image, 
extern unsigned char	BUF_R[LT_3D][X_IM_OP_F][Y_IM_OP_F], //3d input image 
				BUF_G[LT_3D][X_IM_OP_F][Y_IM_OP_F], //in the processing format
				BUF_B[LT_3D][X_IM_OP_F][Y_IM_OP_F];

extern 	int2v im0dms, im1dms,		// storage for dimensions of input images 
		          im2dms;		// aquired from files; currently unused

extern 	bool to_filter;
extern 	double C_extrema_disparity_fltr;

extern 	real4v k1, k2;
extern 	double	PI,		ro_noise, ro_coef,				
		mega1,	mega2,

       shift_per_frm_X, //   shift/frame in simulated sequence
	   shift_per_frm_Y; //

extern 		image_3D     im_sequence;	

extern 		tst_wnd_3D   tw3D;
extern 		th_FM_3D     flow_S, flow, flow_L;
extern 		GLfloat 	colv[4],  
			range_for_visualization,fStep,
			ctrlpoints_R[SCALE_2D][SCALE_2D][3],
			ctrlpoints_G[SCALE_2D][SCALE_2D][3],
			ctrlpoints_B[SCALE_2D][SCALE_2D][3];

extern 		mt19937 eng;  // a core engine cl

void init_im_0( unsigned char * BUFR,
			   unsigned char * BUFG,
			   unsigned char * BUFB, int2v &imdms)
{
	System::Drawing::Bitmap btmp  = Bitmap::FromFile(IM_0_NAME);  /*bad_000.gif gor30.jpg tvety01.jpg kot08.jpg  53_14_sm.tif taxi01.pgm kuhnya8.jpg */
	unsigned char *wrkR, *wrkG, *wrkB; int i,j, bright=0, dark=0, points=0;
	
	double max=-1., min=256./*,  supp_bright, supp_dark*/;
	double wR,wG,wB;
	imdms.y=btmp.Height;  imdms.x=btmp.Width;
	
for( i=0; i<btmp.Width; i++)
		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
			  wrkB=BUFB+i*Y_IM_OP_F; j<btmp.Height; j++)
	{	
		System::Drawing::Color picColor=btmp.GetPixel(i,j);
		wR=(double)picColor.R;
		wR=std_rng_to_mean_to_photons_to_noise(wR);
		if(GREY==true)
		{
		/* grayscale, for example Yosemite*/
		wG=wR; wB=wR;
		}
		else
		/* color, non-Yosemite*/
		{wG=(double)picColor.G;
		wG=std_rng_to_mean_to_photons_to_noise(wG);
		wB=(double)picColor.B;
		wB=std_rng_to_mean_to_photons_to_noise(wB);
		}
		
		*(wrkR+j)= wR+.5 /*picColor.R*/;
		*(wrkG+j)= wG+.5 /*picColor.G*/;
		*(wrkB+j)= wB+.5 /*picColor.B*/;
		Image_0.Im[i][j][0] = (GLubyte)(wR+.5) /*picColor.R*/;
		Image_0.Im[i][j][1] = (GLubyte)(wG+.5) /*picColor.G*/;
		Image_0.Im[i][j][2] = (GLubyte)(wB+.5) /*picColor.B*/;
			//*(wrkR+j)= picColor.R;
			//*(wrkG+j)= picColor.G;
			//*(wrkB+j)= picColor.B;
			////current=0.;
			//Image_0.Im[i][j][0] = (GLubyte)picColor.R;
			////current+=(double)picColor.R*(double)picColor.R;
			//Image_0.Im[i][j][1] = (GLubyte)picColor.G;
			////current+=(double)picColor.G*(double)picColor.G;
			//Image_0.Im[i][j][2] = (GLubyte)picColor.B;
			///* current+=(double)picColor.B*(double)picColor.B;
			//if(max<current)max=current;
			//if(min>current)min=current;
			//if(current>210.*210.)
			//{bright++;}
			//if(current<1.*1.)
			//	dark++; */
			Image_0.Im[i][j][3] = (GLubyte) 255;
	}
	/*supp_bright=(double)bright/(double)(points);
	supp_dark=(double)dark/(double)(points);
	max=sqrt(max);
	min=sqrt(min);*/
	return;
 }

void init_im_1( unsigned char * BUFR,
				unsigned char * BUFG,
				unsigned char * BUFB, int2v &imdms)
{
	System::Drawing::Bitmap btmp  = Bitmap::FromFile(IM_1_NAME); /*bad_001.gif tvety03.jpg kot09.jpg gor31.jpg 53_15_sm.tif taxi02.pgm kuhnya8_25.jpg */
	
	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
	
imdms.y=btmp.Height;
imdms.x=btmp.Width;
double wR,wG,wB;
for( i=0; i<btmp.Width; i++)
		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
			  wrkB=BUFB+i*Y_IM_OP_F; j<btmp.Height; j++)
	{
		System::Drawing::Color picColor=btmp.GetPixel(i,j);
		wR=(double)picColor.R;
		wR=std_rng_to_mean_to_photons_to_noise(wR);

		if(GREY==true)
		{
		/* grayscale, for example Yosemite*/
		wG=wR; wB=wR;
		}
		else
		/* color, non-Yosemite*/
		{wG=(double)picColor.G;
		wG=std_rng_to_mean_to_photons_to_noise(wG);
		wB=(double)picColor.B;
		wB=std_rng_to_mean_to_photons_to_noise(wB);
		}
				
		*(wrkR+j)= wR+.5 /*picColor.R*/;
		*(wrkG+j)= wG+.5 /*picColor.G*/;
		*(wrkB+j)= wB+.5 /*picColor.B*/;
		Image_1.Im[i][j][0] = (GLubyte)(wR+.5) /*picColor.R*/;
		Image_1.Im[i][j][1] = (GLubyte)(wG+.5) /*picColor.G*/;
		Image_1.Im[i][j][2] = (GLubyte)(wB+.5) /*picColor.B*/;
		Image_1.Im[i][j][3] = (GLubyte) 255;
	}
}

void init_im_2( unsigned char * BUFR,
			   unsigned char * BUFG,
			   unsigned char * BUFB, int2v &imdms)
{
	System::Drawing::Bitmap btmp  = Bitmap::FromFile(IM_2_NAME); /*bad_002.gif tvety05.jpg kot10.jpg gor32.jpg 53_16_sm.tif taxi03.pgm kuhnya8_5.jpg */
	double wR,wG,wB;
	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
	
imdms.y=btmp.Height;
imdms.x=btmp.Width;
for( i=0; i<btmp.Width; i++)
		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
			  wrkB=BUFB+i*Y_IM_OP_F; j<btmp.Height; j++)
	{
		System::Drawing::Color picColor=btmp.GetPixel(i,j);
		wR=(double)picColor.R;
		wR=std_rng_to_mean_to_photons_to_noise(wR);
		
		if(GREY==true)
		{
		/* grayscale, for example Yosemite*/
		wG=wR; wB=wR;
		}
		else
		/* color, non-Yosemite*/
		{wG=(double)picColor.G;
		wG=std_rng_to_mean_to_photons_to_noise(wG);
		wB=(double)picColor.B;
		wB=std_rng_to_mean_to_photons_to_noise(wB);
		}
		
		*(wrkR+j)= wR+.5 /*picColor.R*/;
		*(wrkG+j)= wG+.5 /*picColor.G*/;
		*(wrkB+j)= wB+.5 /*picColor.B*/;
		Image_2.Im[i][j][0] = (GLubyte)(wR+.5) /*picColor.R*/;
		Image_2.Im[i][j][1] = (GLubyte)(wG+.5) /*picColor.G*/;
		Image_2.Im[i][j][2] = (GLubyte)(wB+.5) /*picColor.B*/;
		Image_2.Im[i][j][3] = (GLubyte) 255;
	}
}

void init_im_3( unsigned char * BUFR,
			   unsigned char * BUFG,
			   unsigned char * BUFB, int2v &imdms)
{
	System::Drawing::Bitmap btmp  = Bitmap::FromFile(IM_3_NAME); /*bad_003.gif tvety07.jpg kot11.jpg gor33.jpg taxi03.pgm kuhnya8_5.jpg 53_16_sm.tif*/
	double wR,wG,wB;
	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
	
imdms.y=btmp.Height;
imdms.x=btmp.Width;
for( i=0; i<btmp.Width; i++)
		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
			  wrkB=BUFB+i*Y_IM_OP_F; j<btmp.Height; j++)
	{
		System::Drawing::Color picColor=btmp.GetPixel(i,j);
		wR=(double)picColor.R;
		wR=std_rng_to_mean_to_photons_to_noise(wR);

		if(GREY==true)
		{
		/* grayscale, for example Yosemite*/
		wG=wR; wB=wR;
		}
		else
		/* color, non-Yosemite*/
		{wG=(double)picColor.G;
		wG=std_rng_to_mean_to_photons_to_noise(wG);
		wB=(double)picColor.B;
		wB=std_rng_to_mean_to_photons_to_noise(wB);
		}
		
		*(wrkR+j)= wR+.5 /*picColor.R*/;
		*(wrkG+j)= wG+.5 /*picColor.G*/;
		*(wrkB+j)= wB+.5 /*picColor.B*/;
		Image_2.Im[i][j][0] = (GLubyte)(wR+.5) /*picColor.R*/;
		Image_2.Im[i][j][1] = (GLubyte)(wG+.5) /*picColor.G*/;
		Image_2.Im[i][j][2] = (GLubyte)(wB+.5) /*picColor.B*/;
		Image_2.Im[i][j][3] = (GLubyte) 255;
	}
}

//void init_im_0( unsigned char * BUFR,
//			   unsigned char * BUFG,
//			   unsigned char * BUFB, int2v &imdms)
//{
//	/*System::Drawing::Bitmap btmp  = Bitmap::FromFile("yos_img_09.tif"); */ /*bad_000.gif gor30.jpg tvety01.jpg kot08.jpg  53_14_sm.tif taxi01.pgm kuhnya8.jpg */
//	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
//	yos_frame shot("data09");
//	imdms.y=shot.Height;  imdms.x=shot.Width;
//	
//for( i=0; i<shot.Width; i++)
//		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
//			  wrkB=BUFB+i*Y_IM_OP_F; j<shot.Height; j++)
//		{		unsigned char picColor=shot.frame[i][j];
//			*(wrkR+j)= picColor;
//			*(wrkG+j)= picColor;
//			*(wrkB+j)= picColor;
//			Image_0.Im[i][j][0] = (GLubyte)picColor;
//			Image_0.Im[i][j][1] = (GLubyte)picColor;
//			Image_0.Im[i][j][2] = (GLubyte)picColor;
//			Image_0.Im[i][j][3] = (GLubyte) 255;
//	}
// }
//
//void init_im_1( unsigned char * BUFR,
//				unsigned char * BUFG,
//				unsigned char * BUFB, int2v &imdms)
//{
//	/*System::Drawing::Bitmap btmp  = Bitmap::FromFile("yos_img_09.tif"); */ /*bad_000.gif gor30.jpg tvety01.jpg kot08.jpg  53_14_sm.tif taxi01.pgm kuhnya8.jpg */
//	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
//	yos_frame shot("data10");
//	imdms.y=shot.Height;  imdms.x=shot.Width;
//	
//for( i=0; i<shot.Width; i++)
//		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
//			  wrkB=BUFB+i*Y_IM_OP_F; j<shot.Height; j++)
//		{		unsigned char picColor=shot.frame[i][j];
//			*(wrkR+j)= picColor;
//			*(wrkG+j)= picColor;
//			*(wrkB+j)= picColor;
//			Image_1.Im[i][j][0] = (GLubyte)picColor;
//			Image_1.Im[i][j][1] = (GLubyte)picColor;
//			Image_1.Im[i][j][2] = (GLubyte)picColor;
//			Image_1.Im[i][j][3] = (GLubyte) 255;
//	}
//}
//
//void init_im_2( unsigned char * BUFR,
//			   unsigned char * BUFG,
//			   unsigned char * BUFB, int2v &imdms)
//{
//	/*System::Drawing::Bitmap btmp  = Bitmap::FromFile("yos_img_09.tif"); */ /*bad_000.gif gor30.jpg tvety01.jpg kot08.jpg  53_14_sm.tif taxi01.pgm kuhnya8.jpg */
//	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
//	yos_frame shot("data11");
//	imdms.y=shot.Height;  imdms.x=shot.Width;
//	
//for( i=0; i<shot.Width; i++)
//		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
//			  wrkB=BUFB+i*Y_IM_OP_F; j<shot.Height; j++)
//		{		unsigned char picColor=shot.frame[i][j];
//			*(wrkR+j)= picColor;
//			*(wrkG+j)= picColor;
//			*(wrkB+j)= picColor;
//			Image_2.Im[i][j][0] = (GLubyte)picColor;
//			Image_2.Im[i][j][1] = (GLubyte)picColor;
//			Image_2.Im[i][j][2] = (GLubyte)picColor;
//			Image_2.Im[i][j][3] = (GLubyte) 255;
//	}
//	
//}
//
//void init_im_3( unsigned char * BUFR,
//			   unsigned char * BUFG,
//			   unsigned char * BUFB, int2v &imdms)
//{
//	/*System::Drawing::Bitmap btmp  = Bitmap::FromFile("yos_img_09.tif"); */ /*bad_000.gif gor30.jpg tvety01.jpg kot08.jpg  53_14_sm.tif taxi01.pgm kuhnya8.jpg */
//	unsigned char *wrkR, *wrkG, *wrkB; int i,j;
//	yos_frame shot("data12");
//	imdms.y=shot.Height;  imdms.x=shot.Width;
//	
//for( i=0; i<shot.Width; i++)
//		for( j=0, wrkR=BUFR+i*Y_IM_OP_F, wrkG=BUFG+i*Y_IM_OP_F,
//			  wrkB=BUFB+i*Y_IM_OP_F; j<shot.Height; j++)
//		{		unsigned char picColor=shot.frame[i][j];
//			*(wrkR+j)= picColor;
//			*(wrkG+j)= picColor;
//			*(wrkB+j)= picColor;
//			Image_2.Im[i][j][0] = (GLubyte)picColor;
//			Image_2.Im[i][j][1] = (GLubyte)picColor;
//			Image_2.Im[i][j][2] = (GLubyte)picColor;
//			Image_2.Im[i][j][3] = (GLubyte) 255;
//	}
//	
//}


void init_flow_S(image_3D * p_im_sequence, int lt, int ln)
{
	tw3D=    tst_wnd_3D(lt, ln);
	flow_S=  th_FM_3D(tw3D, p_im_sequence);
	
	flow_S.scope=1.;		// it is a dummy parameter for this version
	flow_S.scopeTW=1.;	// it is a dummy parameter for this version
	flow_S.trans_im_3D();  // it is a dummy action for this version
	flow_S.RNGS=ranges(SHIFT_RANGE);       // shift_r_a_n_g_e=10.; 
}
void init_flow(image_3D * p_im_sequence, int lt, int ln)
{
	tw3D=    tst_wnd_3D(lt, ln);
	th_FM_3D temp(tw3D, p_im_sequence);
	flow=   temp;      
	
    flow.scope=1.;		// it is a dummy parameter for this version
	flow.scopeTW=1.;	// it is a dummy parameter for this version
	flow.trans_im_3D();  // it is a dummy action for this version
	flow.RNGS=ranges(SHIFT_RANGE);       // shift_r_a_n_g_e=10.; 
}
void init_flow_L(image_3D * p_im_sequence, int lt, int ln)
{
	tw3D=    tst_wnd_3D(lt, ln);
	flow_L=  th_FM_3D(tw3D, p_im_sequence);
	
	flow_L.scope=1.;		// it is a dummy parameter for this version
	flow_L.scopeTW=1.;	// it is a dummy parameter for this version
	flow_L.trans_im_3D();  // it is a dummy action for this version
	flow_L.RNGS=ranges(SHIFT_RANGE);       // shift_r_a_n_g_e=10.; 
}

void init_var_3D()
{
	poi=simula_Poisson_noise (RO, WCOMMON, BCOMMON); 

	init_im_0( &BUF_R[0][0][0],&BUF_G[0][0][0],&BUF_B[0][0][0], im0dms);
	if(im0dms.x!=X_IM_OP_F||im0dms.y!=Y_IM_OP_F)
		printf("dimension_error\n");
	init_im_1( &BUF_R[1][0][0],&BUF_G[1][0][0],&BUF_B[1][0][0], im0dms);
	if(im0dms.x!=X_IM_OP_F||im0dms.y!=Y_IM_OP_F)
		printf("dimension_error\n");
	init_im_2( &BUF_R[2][0][0],&BUF_G[2][0][0],&BUF_B[2][0][0], im0dms);
	if(im0dms.x!=X_IM_OP_F||im0dms.y!=Y_IM_OP_F)
		printf("dimension_error\n");
	init_im_3( &BUF_R[3][0][0],&BUF_G[3][0][0],&BUF_B[3][0][0], im0dms);
	if(im0dms.x!=X_IM_OP_F||im0dms.y!=Y_IM_OP_F)
		printf("dimension_error\n");

	to_filter=false/*true*/;
	C_extrema_disparity_fltr=C_E_D_th;
	
	double	skewX=0., 
			skewY=0.; 
	int nmb_of_shts =LT_3D, 
		imx=X_IM_OP_F,	imy=Y_IM_OP_F,	imt=LT_3D;

	k1.u=cos(3.*PI/10.);   k1.v=sin(3.*PI/10.);   //  3.*PI/10.=54°
	k2.u=cos(3.*PI/20.);   k2.v=- sin(3.*PI/20.); //  3.*PI/20.=27°
	/*    12mega1/DVA_PI=dx*k1.u + dy*k1.v
	      18mega2/DVA_PI=dx*k2.u + dy*k2.v	

				k2.v*12mega1/DVA_PI=dx*k1.u*k2.v + dy*k1.v*k2.v
				k1.v*18mega2/DVA_PI=dx*k2.u*k1.v + dy*k2.v*k1.v

				k2.v*12mega1/DVA_PI- k1.v*18mega2/DVA_PI= dx*(k1.u*k2.v - k2.u*k1.v)


				k2.u *12mega1/DVA_PI=dx*k1.u*k2.u  + dy*k1.v*k2.u 
				k1.u *18mega2/DVA_PI=dx*k2.u*k1.u +  dy*k2.v*k1.u	

				k2.u *12mega1/DVA_PI - k1.u *18mega2/DVA_PI= dy*(k1.v*k2.u -k2.v*k1.u)
	*/
shift_per_frm_X= (k2.v*12*mega1/DVA_PI- k1.v*18*mega2/DVA_PI)
                                       /((k1.u*k2.v - k2.u*k1.v)*(LT_3D-1));
shift_per_frm_Y= (k2.u*12*mega1/DVA_PI- k1.u*18*mega2/DVA_PI)
                                       /((k1.v*k2.u - k2.v*k1.u)*(LT_3D-1));

/* ro_coef must belong to [0.0005,1], this is  
a noise lelel control for sinusoidal input*/
ro_coef=RO_COEF ;     
//ro_noise=ro_coef*ro_max; 
//Bro_noise=1./(ro_noise*ro_noise);
//Wro_noise=AMD*Bro_noise;

	//poi=simula_Poisson_noise(ro_coef, RO_MX, DMX, DMN); 
	//synthez_im_1_2_3( nmb_of_shts );

	//synthez_3_by_1_shot( shift_per_frm_X, shift_per_frm_Y, skewX, skewY);
	
  	im_sequence= image_3D(imx, imy, imt, &BUF_R[0][0][0],&BUF_G[0][0][0],&BUF_B[0][0][0]);
 
	//AT.init_flow_S(&im_sequence, LT_3D, LN_3DS);
	//AT.init_flow(&im_sequence,   LT_3D, LN_3D);
	//AT.init_flow_L(&im_sequence, LT_3D, LN_3DL);
	

	printf("SCALE: coarse %d  fine %d  super_fine %d \n sequence length %d\n\n",  
		SCALE_2D, SCALE_FINE_GRID_2D, SCALE_SPR_FN_GRID_2D, LT_3D);
	printf("fine_shft_rng=[-%g %g];  2.*shift_r_a_n_g_e/(SCALE_2D-1)=%g\n",  
	   flow.RNGS.vlcty_rng_fine, flow.RNGS.vlcty_rng_fine,
	   2.*flow.RNGS.shift_r_a_n_g_e/(SCALE_2D-1));
	printf("super_fine_shft_rng=[-%g %g]\n\n",  
	   flow.RNGS.rng_spr_fn, flow.RNGS.rng_spr_fn);

    printf("shft_/_frm_X=%g  shft_/_frm_Y=%g  shft_rng=[-%g %g]\n\n noise_level=%g, noise_max_level %g\n", 
					shift_per_frm_X, shift_per_frm_Y, flow.RNGS.shift_r_a_n_g_e, 
					flow.RNGS.shift_r_a_n_g_e, 		poi.ro/*ro_noise*/, poi.ro_max);
	printf("B_ and Wmax=%g  %g  \n sigma_B_ and Wmax=%g %g\n", poi.Bmax, poi.Wmax,
					poi.sigma_Bmax, poi.sigma_Wmax);
	printf("B_ and Wcommon=%g  %g  \n B_ and Wro_noise=%g %g\n\n", poi.Bcommon, poi.Wcommon,
					poi.Bro_noise, poi.Wro_noise);


//test on coordinated creation of flow.IM (for image processing)
//and Image_* (for visualization of input)

//trns_BUF_to_Image(); 

//int m, n;
//     // Initialize spline variable
//
//	m=SCALE_2D;
//	n=SCALE_2D;
   
     GLint iLine, jLine;
	int i,j;
	/* while scanning in 'scan_OF_cand' 
	del=(2.*range/(float)(scale_2D-1)); hence, range=(scale_2D-1)del/2   */
	
	range_for_visualization=0.5f*fStep*(SCALE_2D-1); 
	  
	colv[3]=0.75f;
		
	for(i=0, iLine= -range_for_visualization; i<SCALE_2D; i++, iLine += fStep)
		for(j=0, jLine= -range_for_visualization; j<SCALE_2D;  j++, jLine += fStep)
		{	ctrlpoints_R[i][j][0]=iLine;       /* x-coordinate  of plot */
			ctrlpoints_R[i][j][2]=jLine;       /* z-coordinate   */

			ctrlpoints_G[i][j][0]=iLine;	
			ctrlpoints_G[i][j][2]=jLine;

	        ctrlpoints_B[i][j][0]=iLine;	
			ctrlpoints_B[i][j][2]=jLine;	
		}

}

void synthez_3_by_1_shot(double shift_per_frm_X, double shift_per_frm_Y, double skewX, double skewY)
{

	init_im_1( &BUF_R[1][0][0],&BUF_G[1][0][0],&BUF_B[1][0][0], im1dms);

		int i, j, time; 

	for (time=0; time<3; time+=2)
	for (i=0; i<X_IM_OP_F; i++)
		for (j=0; 	j<Y_IM_OP_F; j++)
		{	
			double d=(double)(LT_3D-1), 
			dtime=((double)time - d*0.5)/d;

			int i0= i+ dtime*shift_per_frm_X +0.5+ ABS(dtime)*skewX,
				j0= j+ dtime*shift_per_frm_Y +0.5+ ABS(dtime)*skewY;

			if(i0<0)i0=0; if(i0>=X_IM_OP_F)i0=X_IM_OP_F-1;
			if(j0<0)j0=0; if(j0>=Y_IM_OP_F)j0=Y_IM_OP_F-1;

			
					BUF_R[time][i0][j0]=BUF_R[1][i][j];
					BUF_G[time][i0][j0]=BUF_G[1][i][j];
					BUF_B[time][i0][j0]=BUF_B[1][i][j];
			
				switch(time){
						case 0:{
					Image_0.Im[i0][j0][0] = (GLubyte)BUF_R[1][i][j];
					Image_0.Im[i0][j0][1] = (GLubyte)BUF_G[1][i][j];
					Image_0.Im[i0][j0][2] = (GLubyte)BUF_B[1][i][j];
					Image_0.Im[i0][j0][3] = (GLubyte) 255; break;
							   }
					   case 2:{
					Image_2.Im[i0][j0][0] = (GLubyte)BUF_R[1][i][j];
					Image_2.Im[i0][j0][1] = (GLubyte)BUF_G[1][i][j];
					Image_2.Im[i0][j0][2] = (GLubyte)BUF_B[1][i][j];
					Image_2.Im[i0][j0][3] = (GLubyte) 255; break;
							   }
				}
}

}

double brightness(int i, int j, int time)
{ 
	double fi,fj,  r, mean;
		//mega1=3./PI, mega2=3./DVA_PI;

	fi=i;	fj=j;

	/*real4v k1, k2;

	k1.u=cos(3.*PI/10.);   k1.v=sin(3.*PI/10.); 
	k2.u=cos(3.*PI/20.);   k2.v=- sin(3.*PI/20.); */

	 //noise level  is presented by 'ro_noise'; it belongs to [0; ro_max]
      	
	r=0.25*(2.+   sin(DVA_PI*(k1.u*fi+k1.v*fj)/12.  +mega1*time/*/(LT_3D-1)*/)+ 
				  sin(DVA_PI*(k2.u*fi+k2.v*fj)/18.  +mega2*time/*/(LT_3D-1)*/) );
    
	//mean=Bro_noise+(Wro_noise-Bro_noise)*r;
	mean=poi.EXPECT_Photons(r);
    poisson_distribution<double> poisson(mean);
	//r= (1.-2.*ro_noise)*r+ro_noise;             //  ro_noise<r<1-ro_noise
	//noise= ro_noise*(2.*(double)rand() / (double)(RAND_MAX + 1) - 1.); // -ro_noise<noise<ro_noise

	return poisson(eng);
	
};

void synthez_im_1_2_3(  int nmb_of_shts )
{   
	/*  !!! this syntesis is valid for greyscle only, 
	    be careful with repsect to color noise-components 
		of 'brightness' !!! */

int i,j, time;
unsigned char *wrkR, *wrkG, *wrkB;
for (time=0; time<nmb_of_shts; time++)
	for (i=0; i<X_IM_OP_F; i++)
		for (j=0, wrkR=&BUF_R[time][0][0]/*BUFR*/+i*Y_IM_OP_F,
				  wrkG=&BUF_G[time][0][0]+i*Y_IM_OP_F,
				  wrkB=&BUF_B[time][0][0]+i*Y_IM_OP_F; 
										j<Y_IM_OP_F; j++)
		{        double br=  brightness(i,j,time);

			// begin: corrected signal to fit Poisson model
				br=poi.map_Photons_to_std_range(br);
				/*br=Bcommon+(Wcommon-Bcommon)*(br-Bro_noise)/(Wro_noise-Bro_noise);
				if(br<0.)br=0.;
				if(br>255.)br=255.;*/
				br/=255;	// this operation is introduced to avoid corrections 
							// in the following texts of a previous version
			// end: corrected signal to fit Poisson model

				*(wrkR+j)= (unsigned char) (br* 255. +0.5);
				*(wrkG+j)= (unsigned char) (br* 255.+0.5);
				*(wrkB+j)= (unsigned char) (br* 255.+0.5);
				switch(time){
						case 0:{
					Image_0.Im[i][j][0] = (GLubyte)(br* 255.+0.5);
					Image_0.Im[i][j][1] = (GLubyte)(br* 255.+0.5);
					Image_0.Im[i][j][2] = (GLubyte)(br* 255.+0.5);
					Image_0.Im[i][j][3] = (GLubyte) 255; break;
							   }
						case 1:{
					Image_1.Im[i][j][0] = (GLubyte)(br* 255.+0.5);
					Image_1.Im[i][j][1] = (GLubyte)(br* 255.+0.5);
					Image_1.Im[i][j][2] = (GLubyte)(br* 255.+0.5);
					Image_1.Im[i][j][3] = (GLubyte) 255; break;
							   }
					   case 2:{
					Image_2.Im[i][j][0] = (GLubyte)(br* 255.+0.5);
					Image_2.Im[i][j][1] = (GLubyte)(br* 255.+0.5);
					Image_2.Im[i][j][2] = (GLubyte)(br* 255.+0.5);
					Image_2.Im[i][j][3] = (GLubyte) 255; break;
							   }
					   default: break;
				}
}
}