/******************************************************************************************
> [KAZOE]                                                                                 <
> Initial Sofware by [Yannick Verdie, Florent Lafarge], Copyright (C) Inria (GEOMETRICA)  <
> All Rights Reserved, 2012, [Version 1.0].                                               <
******************************************************************************************/
#include "GPU_Flamengo.cuh"
#include "DataTree.h"



__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Death_on_Grid_RJMCMC(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							  const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x ;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];


	if (idx <  nb_Cells && index_obj >= 0)
	{
		//--------------------------------get ellipse----------------------------------
		const float CenterX		= result_Ellipses->dev_x[index_obj];
		const float CenterY		= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 
		//---------------------------------------------------------------------------------

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;									//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;									//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		float resInter = 0;
		unsigned int nbNewPixels_death = 0;float sigmaEbackground_death = 0; float sigmaEcircle_death = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed

		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;


		for (int y = -searchArea ; y <= searchArea; y++ )
		{
			for (int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f );
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f );

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f ) area++;
			}
		}

		for (int  y = -searchArea ; y <= searchArea; y++ )
		{	  
			for (int  x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f );
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f );


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f && coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) //in ellipse and in image with border
					{

						tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// and in image
						{
							current_val = val[tmp1];
							if (nb_intersection == 1)//death
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg);
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}

						if ( nb_intersection > 1)//	(1 is itself)
						{
							bool already_jump = false;
							tmp1 = area;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area && !already_jump)
								{
									already_jump = true;
								}else{
									tmp1 = max(tmp1,val_area);
								}
							}
							resInter += (1.f / min ( area , tmp1));
						}

					}
				}
			}
		}

		if ( resInter > param->Hard_Core ) resInter = 1.0f;

		float sigma_death = -resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha * 
			(
			+ sigmaEbackground_death / (2.f * param->Background_Sigma) - nbNewPixels_death * 0.5f * logf   ( 2.f * CUDART_PI_F * param->Background_Sigma )
			-     sigmaEcircle_death / (2.f * param->Shape_Sigma     ) + nbNewPixels_death * 0.5f * logf   ( 2.f * CUDART_PI_F * param->Shape_Sigma      ) 
			);

	   	tmp1 = global_x_origin + 2*nb_cell_w*global_y_origin;
		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[tmp1];


		if ( evaluateGreen(sigma_death  , temperature ,nb_Ellipse_in_this_cell/float(local_size_cell_w*local_size_cell_h), dev_rando[idx]) )
		{
			//delete ellipse
			result_Ellipses->dev_available_spot[index_obj] = true;
			dev_nbObjectsAtCell_i[tmp1]	= nb_Ellipse_in_this_cell - 1;

			for (int y = -searchArea ; y <= searchArea; y++ )
			{
				for (int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = ((axis0vec[0]*x)*(axis0vec[0]*x)+(axis0vec[1]*y)*(axis0vec[1]*y)) / (axis0 + 0.5f );
					ratio1 = ((axis0vec[1]*x)*(axis0vec[1]*x)+(-axis0vec[0]*y)*(-axis0vec[0]*y)) / (axis1 + 0.5f );

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);
	
						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) 
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];
	
							//if (nb_intersection > 0) //not possible 'false'
							{
								for (int index2 = 0;index2 < nb_intersection-1 ; index2++)
								{	
									if (dev_minArea[tmp1+(w*h)*index2] == area)
									{
										dev_minArea[tmp1+(w*h)*index2] = dev_minArea[tmp1+(w*h)*(nb_intersection-1)];//swap last one with this one
										break;
									}
								}

								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection-1;
							}
						}

					}
				}
			} 
		}
	}
}


__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Death_on_Grid_RJMCMC_bhat(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, const float ratio_birth_death,								//output
							   const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{
	
	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell

	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];

	if (idx <  nb_Cells && index_obj >= 0)
	{
		//--------------------------------get ellipse----------------------------------
		const float CenterX		= result_Ellipses->dev_x[index_obj];
		const float CenterY		= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 
		//---------------------------------------------------------------------------------

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;									//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;									//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell


		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0;

		float resInter = 0;
		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;
		const int border_boundary = param->Bhatt_Thickness;

		const int searchAreaExt = searchArea + border_boundary;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;


		for (int y = -searchArea ; y <= searchArea; y++ )
		{
			for (int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int tmp1, tmp2;
		unsigned int val_area;
		int current_val;//need to be signed because we substract it to something else
		unsigned int nb_intersection;
		bool in_small, in_big, in_space, in_image;
		bool already_jump;
		for (int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{
			for (int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

					if (in_small && in_space)//in ellipse and image with border
					{

						tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// in image
						{
							current_val = val[tmp1];
							nbNewPixelsm1++;
							sigmaEcirclem1      += current_val * current_val;
							meanEcirclem1		+= current_val;
						}

						if ( nb_intersection > 1)//	(1 is itself)
						{
							already_jump = false;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area && !already_jump)
								{
									already_jump = true;
								}else{					
									resInter += (1.0f/ min ( area , val_area));
								}
							}
						}							
					}

					if (in_big && in_image && !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		float sigma_death = 0;

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core ) resInter = 1.0f;

			sigma_death = -getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);		
		}

		tmp1 = global_x_origin + 2*nb_cell_w*global_y_origin;
		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[tmp1];

		if (resInter > param->Hard_Core || (nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen(sigma_death, temperature ,(ratio_birth_death/(1-ratio_birth_death))*nb_Ellipse_in_this_cell/float(local_size_cell_w*local_size_cell_h), dev_rando[idx])) )
		{
			//delete ellipse
			result_Ellipses->dev_available_spot[index_obj] = true;
			dev_nbObjectsAtCell_i[tmp1] = nb_Ellipse_in_this_cell - 1;


			for (int y = -searchArea ; y <= searchArea; y++ )
			{
				for (int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);
	
						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) 
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							//if (nb_intersection > 0)		//not possible
							{
								for (int index2 = 0;index2 < nb_intersection-1 ; index2++)	  //the last one is delete
								{	
									val_area = dev_minArea[tmp1+(w*h)*index2];
									if (val_area == area)
									{
										dev_minArea[tmp1+(w*h)*index2] = dev_minArea[tmp1+(w*h)*(nb_intersection-1)];//swap last one with this one
										break;
									}
								}

								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection-1;
							}
						}

					}
				}
			} 
		}
	}
}




//new	 but not used, slower than the implementation bellow (called old)
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Translate_on_Grid2(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x) ;//idx is the index of the first pixel of the grill cell

	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];


	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 

		float newCenterX; float newCenterY;
		unsigned int area = 0;

		//1sd method (uniform ?)
		float theta = dev_rando[idx] * 2 * CUDART_PI_F;
		float rho =	 dev_rando[idx + nb_Cells] * param -> MaxTranslation;
		newCenterX =CenterX + cosf(theta)*rho;
		newCenterY =CenterY + sinf(theta)*rho ;

		//stay in space
		if (newCenterX < -border || newCenterX >= w - border || newCenterY < -border || newCenterY >= h - border) return;

		float sigma = 0.0f;
		  //Cost to kill the current shape
		 float resInter = 0;unsigned int nbNewPixels_death = 0;float sigmaEbackground_death = 0; float sigmaEcircle_death = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		area = 0;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);

		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int position_area;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{	  
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{

				//cost kill + remove
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) //in ellipse and in image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 1)//death
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg);
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}


						if ( nb_intersection > 1)//	(1 is itself)
						{
							 position_area = -1;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area && position_area < 0)
								{
									position_area = index2;
								}else{
									resInter += (1.0f/ min ( area , val_area));
								}
							}
						}
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;

		sigma = -resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha * 
			(
			+ sigmaEbackground_death / (2.f* param->Background_Sigma) - nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Background_Sigma )
			-     sigmaEcircle_death / (2.f* param->Shape_Sigma     ) + nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Shape_Sigma      ) 
			);
		//-------------------------------------------end cost to kill previous


		bool in_both;bool isInside = false;
		resInter = 0;nbNewPixels_death = 0;sigmaEbackground_death = 0; sigmaEcircle_death = 0;
		//-------------------------------------------cost to create new
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = newCenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = newCenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
					{					
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !

						ratio0 = (axis0vec[0]*(coordX-CenterX) + axis0vec[1]*(coordY-CenterY)) / (axis0 + 0.5f);
						ratio1 = (axis0vec[1]*(coordX-CenterX) - axis0vec[0]*(coordY-CenterY)) / (axis1 + 0.5f);
						in_both = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];
						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// in image
						{
							isInside = true;
							tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
							current_val = val[tmp1];
							if (nb_intersection <= in_both)	//0 if the new does not intersect the old one (and any other) or 1 if the new does intersect the old one (only)
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg) ;
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}
						
					
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area , val_area));
						}

						if (in_both)  resInter -= (1.0f/area);
						
						
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;


		sigma += resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha *
			(
			- sigmaEbackground_death / (2.f* param->Background_Sigma) + nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Background_Sigma )
			+     sigmaEcircle_death / (2.f* param->Shape_Sigma )     - nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Shape_Sigma      )
			);

	
	   if (isInside && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+2*nb_Cells]))
       {
		   //unsigned	int nb_overlap;
		   unsigned int val_area;
			
			//remove previous
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);


 
						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)
						{
						    tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];
							
							//if ( nb_intersection > 0)
							{
								for (int index2 = 0;index2 < nb_intersection-1 ; index2++)
								{	
									val_area = dev_minArea[tmp1+(w*h)*index2];
									if (val_area == area)
									{
										dev_minArea[tmp1+(w*h)*index2] = dev_minArea[tmp1+(w*h)*(nb_intersection-1)];//swap last one with this one
										break;
									}
								}

								
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection-1;		
							}  
							
						}

						  
						//add new one
						coordX = newCenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = newCenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);



						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}   
	
			//update X and Y
			result_Ellipses->dev_x[index_obj] = newCenterX;
			result_Ellipses->dev_y[index_obj] = newCenterY;

			unsigned int global_x_origin = int((CenterX+border)/size_cell);
			unsigned int global_y_origin = int((CenterY+border)/size_cell);
			//!update nb_Ellipse_per_cell and idCell if the object change of cell
			int new_x_id = int((newCenterX+border)/size_cell);
			int new_y_id = int((newCenterY+border)/size_cell);										//index based on grid reference


			if (new_x_id != global_x_origin || new_y_id != global_y_origin )
			{		
				dev_nbObjectsAtCell_i[global_x_origin + int(w/size_cell + 1)*global_y_origin]--;
				dev_nbObjectsAtCell_i[new_x_id + int(w/size_cell + 1)*new_y_id]++;
			}
		}		
	}
}


//kill the old and add it back after   -> BUT FASTER, WHY ?
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Translate_on_Grid(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell	
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];


	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX		= result_Ellipses->dev_x[index_obj];
		const float CenterY		= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;									//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;									//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		float newCenterX; float newCenterY;
		unsigned int area = 0;

		//1sd method (uniform ?)
		float theta = dev_rando[idx] * 2 * CUDART_PI_F;
		float rho =	 dev_rando[idx + nb_Cells] * param -> MaxTranslation;
		newCenterX =CenterX + cosf(theta)*rho;
		newCenterY =CenterY + sinf(theta)*rho ;

		//stay in space
		if (newCenterX < -border || newCenterX >= w - border || newCenterY < -border || newCenterY >= h - border) return;


		float sigma = 0.0f;
		  //Cost to kill the current shape
		 float resInter = 0;unsigned int nbNewPixels_death = 0;float sigmaEbackground_death = 0; float sigmaEcircle_death = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		area = 0;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int position_area;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{	  
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{

				//cost kill + remove
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) //in ellipse and in image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 1)//death
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg);
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}


						if ( nb_intersection > 1)//	(1 is itself)
						{
							 position_area = -1;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area && position_area < 0)
								{
									position_area = index2;
								}else{
									resInter += (1.0f/ min ( area , val_area));
								}
							}
						}

						//remove at the same time
						if (nb_intersection > 0)
						{
							if (nb_intersection > 1)
								dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

							dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
						}

					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;

		sigma = -resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha * 
			(
			+ sigmaEbackground_death / (2.f* param->Background_Sigma) - nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Background_Sigma )
			-     sigmaEcircle_death / (2.f* param->Shape_Sigma     ) + nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Shape_Sigma      ) 
			);
		//-------------------------------------------end cost to kill previous


		bool isInside = false;
		resInter = 0;nbNewPixels_death = 0;sigmaEbackground_death = 0; sigmaEcircle_death = 0;
		//-------------------------------------------cost to create new
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = newCenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = newCenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
					{					
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !

						ratio0 = (axis0vec[0]*(x-CenterX+newCenterX) + axis0vec[1]*(y-CenterY+newCenterY)) / (axis0 + 0.5f);
						ratio1 = (axis0vec[1]*(x-CenterX+newCenterX) - axis0vec[0]*(y-CenterY+newCenterY)) / (axis1 + 0.5f);

						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];
						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// in image
						{
							isInside = true;
							tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
							current_val = val[tmp1];
							if (nb_intersection == 0)	//0 if the new does not intersect the old one (and any other) or 1 if the new does intersect the old one (only)
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg) ;
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}
						
					
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area , val_area));
						} 
						
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;


		sigma += resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha *
			(
			- sigmaEbackground_death / (2.f* param->Background_Sigma) + nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Background_Sigma )
			+     sigmaEcircle_death / (2.f* param->Shape_Sigma )     - nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Shape_Sigma      )
			);

	
	   if (isInside && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+2*nb_Cells]))
       {
			  
			//add new one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = newCenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = newCenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);



						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}   
	
			//update X and Y
			result_Ellipses->dev_x[index_obj] = newCenterX;
			result_Ellipses->dev_y[index_obj] = newCenterY;

			//!update nb_Ellipse_per_cell and idCell if the object change of cell
			int new_x_id = int((newCenterX+border)/size_cell);
			int new_y_id = int((newCenterY+border)/size_cell);	//index based on grid reference


			if (new_x_id != global_x_origin || new_y_id != global_y_origin )
			{	
				dev_nbObjectsAtCell_i[global_x_origin + 2*nb_cell_w*global_y_origin]--;
				dev_nbObjectsAtCell_i[new_x_id + 2*nb_cell_w*new_y_id]++;
			}
		}
		

		else{
			
	 		//add old one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}
		}		
	}
}


   
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Translate_on_Grid_bhat(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];

	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 
		//---------------------------------------------------------------------------------
		

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;	    				//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;	//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		if ( global_x_origin*size_cell >= w || global_y_origin*size_cell >= h) return;

		float newCenterX; float newCenterY;
		unsigned int area = 0;

		float theta = dev_rando[idx] * CUDART_PI_F;
		float rho =	 dev_rando[idx + nb_Cells] * param -> MaxTranslation;
		newCenterX =CenterX + cosf(theta)*rho;
		newCenterY =CenterY + sinf(theta)*rho ;

		//stay in space
		if (newCenterX < -border || newCenterX >= w - border || newCenterY < -border || newCenterY >= h - border) return;


		  //Cost to kill the current shape
		 float resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0; 
		 
		 float sigma = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		const int border_boundary = param->Bhatt_Thickness;
		const int searchAreaExt = searchArea + border_boundary;

		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);

		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		bool in_small, in_big, in_space, in_image;
		int position_area;
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{	  
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{

				//cost kill + remove
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

					if (in_small && in_space)//in ellipse and image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 1)//death
							{
								nbNewPixelsm1++;
								sigmaEcirclem1      += current_val * current_val;
								meanEcirclem1		+= current_val;
							}
						}


						if ( nb_intersection > 1)//	(1 is itself)
						{
							 position_area = -1;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area && position_area < 0)
								{
									position_area = index2;
								}else{
									resInter += (1.0f/ min ( area , val_area));
								}
							}
						}

						//remove at the same time
						if (nb_intersection > 0)
						{
							if (nb_intersection > 1)
								dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

							dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
						}

					}

					if (in_big && in_image && !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma = -getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}
		//-------------------------------------------end cost to kill previous

		resInter = 0;
		nbNewPixelsm1 = 0;meanEcirclem1 = 0; sigmaEcirclem1 = 0;
		nbNewPixelsm2 = 0;meanEcirclem2 = 0; sigmaEcirclem2 = 0; 
		//-------------------------------------------cost to create new
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				{
					coordX = newCenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = newCenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

					if (in_small && in_space)//in ellipse and image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 0)
							{
								nbNewPixelsm1++;
								sigmaEcirclem1      += current_val * current_val;
								meanEcirclem1		+= current_val;
							}
						}

						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area , val_area));
						} 

					}

					if (in_big && in_image && !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma += getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}

		int new_x_id = int((newCenterX+border)/size_cell);
		int new_y_id = int((newCenterY+border)/size_cell);	//index based on grid reference

		val_area = 1;
		if (new_x_id != global_x_origin || new_y_id != global_y_origin )	//change the entensity of the measure if we change the cell, otherwise 1
			val_area = dev_nbObjectsAtCell_i[global_x_origin + 2*nb_cell_w*global_y_origin]/(dev_nbObjectsAtCell_i[ new_x_id + 2*nb_cell_w*new_y_id]+1);
	
	   if (nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma, temperature , val_area , dev_rando[idx+2*nb_Cells]))
       {

			//add new one
			for ( int y = -searchArea; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = newCenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = newCenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);



						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}   
	
			//update X and Y
			result_Ellipses->dev_x[index_obj] = newCenterX;
			result_Ellipses->dev_y[index_obj] = newCenterY;

			//!update nb_Ellipse_per_cell and idCell if the object change of cell
			if (new_x_id != global_x_origin || new_y_id != global_y_origin )
			{
				dev_nbObjectsAtCell_i[global_x_origin + 2*nb_cell_w*global_y_origin]--;
				dev_nbObjectsAtCell_i[new_x_id + 2*nb_cell_w*new_y_id]++;
			}
		}
		

		else{
			
	 		//add old one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}
		}		
	}
}













__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Rotate_on_Grid(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x) ;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];


	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;	    				//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;	//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		if ( global_x_origin*size_cell >= w || global_y_origin*size_cell >= h) return;

		//1sd method (uniform ?)
		float theta = dev_rando[idx] * 2 * CUDART_PI_F;

		float sigma = 0.0f;
		  //Cost to kill the current shape
		 float resInter = 0;
		unsigned int nbNewPixels_death = 0;float sigmaEbackground_death = 0; float sigmaEcircle_death = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		unsigned int area_old = 0;	unsigned int area_new = 0;
		float axis0vec_old[2];axis0vec_old[0]=cosf(angle);axis0vec_old[1]=sinf(angle);
		float axis1vec_old[2];axis1vec_old[0]=axis0vec_old[1];axis1vec_old[1]=-axis0vec_old[0];
		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec_old[0]*x + axis0vec_old[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis1vec_old[0]*x + axis1vec_old[1]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_old++;
			}
		}

		//cost death
		int position_area;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{	  
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec_old[0]*x + axis0vec_old[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis1vec_old[0]*x + axis1vec_old[1]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) //in ellipse and in image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 1)//death
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg);
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}

						if ( nb_intersection > 1)//	(1 is itself)
						{
							position_area = -1;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area_old && position_area < 0)//!already_jump)
								{
									position_area = index2;
								}else{
									resInter += (1.0f/ min ( area_old , val_area));
								}
							}
						}

						//remove at the same time
						if (nb_intersection > 0)
						{
							if (nb_intersection > 1)
								dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

							dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
						}
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;

		sigma = -resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha * 
			(
			+ sigmaEbackground_death / (2.f* param->Background_Sigma) - nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Background_Sigma )
			-     sigmaEcircle_death / (2.f* param->Shape_Sigma     ) + nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Shape_Sigma      ) 
			);

		//-------------------------------------------end cost to kill previous


	    bool isInside = false;
		resInter = 0;nbNewPixels_death = 0;sigmaEbackground_death = 0; sigmaEcircle_death = 0;
		float axis0vec_new[2];axis0vec_new[0]=cosf(theta);axis0vec_new[1]=sinf(theta);
		float axis1vec_new[2];axis1vec_new[0]=axis0vec_new[1];axis1vec_new[1]=-axis0vec_new[0];
		//-------------------------------------------cost to create new
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis1vec_new[0]*x + axis1vec_new[1]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_new++;
			}
		}
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis1vec_new[0]*x + axis1vec_new[1]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
					{
						
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];
						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// in image
						{
							isInside = true;
							if (nb_intersection == 0)	//0 if the new does not intersect the old one (and any other) or 1 if the new does intersect the old one (only)
							{
								tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
								current_val = val[tmp1];
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg) ;
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}
						
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area_new , val_area));
						}				
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;


		sigma += resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha *
			(
			- sigmaEbackground_death / (2.f* param->Background_Sigma) + nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Background_Sigma )
			+     sigmaEcircle_death / (2.f* param->Shape_Sigma )     - nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Shape_Sigma      )
			);

	
	   if (isInside && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+nb_Cells]))
       {
			//add new one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis1vec_new[0]*x + axis1vec_new[1]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);


						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)	//during the remove step, I already update the cell that are common to both ellipse
							{
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_new;
							}
						}
					}
				}
			}   		
			
			//update angle
			result_Ellipses->dev_angle[index_obj] = theta;	 


	   }else{

		 //add old one back
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec_old[0]*x + axis0vec_old[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis1vec_old[0]*x + axis1vec_old[1]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							 nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if ( nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)* nb_intersection] = area_old;
								dev_nbObjectsAtPixel_i[tmp1] =  nb_intersection+1;
							}
						}
						
						
					}
				}
			}
	   }
	}
}





__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Rotate_on_Grid_bhat(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{

	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];

	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 
		//---------------------------------------------------------------------------------
		

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;	    				//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;	//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		if ( global_x_origin*size_cell >= w || global_y_origin*size_cell >= h) return;



		//1sd method (uniform ?)
		float theta;
		if (ELLIPSE_STEP_ANGLE > 0)
			theta   = int(dev_rando[idx] * ELLIPSE_STEP_ANGLE) * CUDART_PI_F / ELLIPSE_STEP_ANGLE ;
		else
			theta   = dev_rando[idx] * CUDART_PI_F ;

		float axis0vec_new[2];axis0vec_new[0]=cosf(theta);axis0vec_new[1]=sinf(theta);

		float sigma = 0.0f;

		//Cost to kill the current shape
		 unsigned int area_old = 0;unsigned int area_new = 0;
		 float resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0; 
		 
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		const int border_boundary = param->Bhatt_Thickness;
		const int searchAreaExt = searchArea + border_boundary;

		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_old++;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_new++;
			}
		}

		bool in_small, in_big, in_space, in_image;
		int position_area;
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{	  
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				//cost kill + remove
				coordX = CenterX + x;
				if ( coordX < 0 )  coordX = (int)(coordX - 1);
				else               coordX = (int)(coordX + 1);

				coordY = CenterY + y;
				if ( coordY < 0 )  coordY = (int)(coordY - 1);
				else               coordY = (int)(coordY + 1);

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
				in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
				in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
				in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

				if (in_small && in_space)//in ellipse and image with border
				{
					tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
					nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

					if ( in_image)// and in image
					{
						tmp1 = coordX+(w-2*border)*(coordY);
						current_val = val[tmp1];//origin at (0,0) image
						nbNewPixelsm1++;
						sigmaEcirclem1      += current_val * current_val;
						meanEcirclem1		+= current_val;
					}


					if ( nb_intersection > 1)//	(1 is itself)
					{
						 position_area = -1;
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							if (val_area == area_old && position_area < 0)
							{
								position_area = index2;
							}else{
								resInter += (1.0f/ min ( area_old , val_area));
							}
						}
					}

					//remove at the same time
					if (nb_intersection > 0 )
					{
						if (nb_intersection > 1 && position_area >= 0)
							dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

						dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
					}

				}

				if (in_big && in_image && !in_small)//in ellipse and image with border
				{
					tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
					current_val = val[tmp1];
					nbNewPixelsm2++;
					sigmaEcirclem2      += current_val * current_val;
					meanEcirclem2		+= current_val;						
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma = -getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}
		//-------------------------------------------end cost to kill previous

		resInter = 0;
		nbNewPixelsm1 = 0;meanEcirclem1 = 0; sigmaEcirclem1 = 0;
		nbNewPixelsm2 = 0;meanEcirclem2 = 0; sigmaEcirclem2 = 0; 
		//-------------------------------------------cost to create new
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				coordX = CenterX + x;
				if ( coordX < 0 )  coordX = (int)(coordX - 1);
				else               coordX = (int)(coordX + 1);

				coordY = CenterY + y;
				if ( coordY < 0 )  coordY = (int)(coordY - 1);
				else               coordY = (int)(coordY + 1);

				ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec_new[1]*x - axis0vec_new[0]*y) / (axis1 + 0.5f);
				in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f+border_boundary);
				ratio1 = (axis0vec_new[1]*x - axis0vec_new[0]*y) / (axis1 + 0.5f+border_boundary);
				in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
				in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

				if (in_small && in_space)//in ellipse and image with border
				{
					tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
					nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

					if ( in_image)// and in image
					{
						tmp1 = coordX+(w-2*border)*(coordY);
						current_val = val[tmp1];//origin at (0,0) image
						nbNewPixelsm1++;
						sigmaEcirclem1      += current_val * current_val;
						meanEcirclem1		+= current_val;
					}

					for (int index2 = 0;index2 < nb_intersection ; index2++)
					{	
						val_area = dev_minArea[tmp2+(w*h)*index2];
						resInter += (1.0f/ min ( area_new , val_area));
					} 

				}

				if (in_big && in_image && !in_small)//in ellipse and image with border
				{
					tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
					current_val = val[tmp1];
					nbNewPixelsm2++;
					sigmaEcirclem2      += current_val * current_val;
					meanEcirclem2		+= current_val;						
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma += getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}

	
	   if (nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+2*nb_Cells]))
       {
			//add new one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec_new[0]*x + axis0vec_new[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec_new[1]*x - axis0vec_new[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);



						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_new;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}   
	
			result_Ellipses->dev_angle[index_obj] = theta; 
		}
		

		else{
			
	 		//add old one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_old;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}
		}		
	}
}






__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void EnlargeShrink_on_Grid(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{
 
	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x) ;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];


	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;	    				//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;	//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		if ( global_x_origin*size_cell >= w || global_y_origin*size_cell >= h) return;



		float newAxis0; float newAxis1;

		//1sd method (uniform ?)
		newAxis0 = param->Ellipse_MinWidth  + dev_rando[idx]				  * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth); 
		newAxis1 = param->Ellipse_MinHeight + dev_rando[idx+ nb_Cells]        * ( param->Ellipse_MaxHeight - param->Ellipse_MinHeight);
		float sigma = newAxis0;//use as temporary
		newAxis0 = max(sigma, newAxis1); //axis 0 is bigger than axis 1
		newAxis1 = min(sigma, newAxis1);
		sigma = 0;

		  //Cost to kill the current shape
		 float resInter = 0;
		unsigned int nbNewPixels_death = 0;float sigmaEbackground_death = 0; float sigmaEcircle_death = 0;
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		unsigned int area_old = 0;	unsigned int area_new = 0;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;//can be signed
		int tmp2;unsigned int val_area;
		int current_val;//must be signed
		unsigned int nb_intersection;

		//get area
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_old++;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_new++;
			}
		}

		int position_area;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{	  
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border) //in ellipse and in image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// and in image
						{
							tmp1 = coordX+(w-2*border)*(coordY);
							current_val = val[tmp1];//origin at (0,0) image
							if (nb_intersection == 1)//death
							{
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg);
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}

						if ( nb_intersection > 1)//	(1 is itself)
						{
							position_area = -1;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = dev_minArea[tmp2+(w*h)*index2];
								if (val_area == area_old && position_area < 0)
								{
									position_area = index2;
								}else{
									resInter += (1.0f/ min ( area_old , val_area));
								}
							}
						}

						//remove at the same time
						if (nb_intersection > 0)
						{
							if (nb_intersection > 1)
								dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

							dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
						}
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 1.0f;

		sigma = -resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha * 
			(
			+ sigmaEbackground_death / (2.f* param->Background_Sigma) - nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Background_Sigma )
			-     sigmaEcircle_death / (2.f* param->Shape_Sigma     ) + nbNewPixels_death * 0.5f* log   ( 2.f* CUDART_PI_F * param->Shape_Sigma      ) 
			);

		//-------------------------------------------end cost to kill previous


	    bool isInside = false;
		resInter = 0;nbNewPixels_death = 0;sigmaEbackground_death = 0; sigmaEcircle_death = 0;
		//-------------------------------------------cost to create new
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = CenterX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = CenterY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
					{	
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];
						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// in image
						{
							isInside = true;
							if (nb_intersection == 0)	//0 if the new does not intersect the old one (and any other) or 1 if the new does intersect the old one (only)
							{
								tmp1 = coordX+(w-2*border)*(coordY);//origin at (0,0) image
								current_val = val[tmp1];
								nbNewPixels_death++;
								sigmaEbackground_death  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg) ;
								sigmaEcircle_death      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}
						
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area_new , val_area));
						}				
					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 10000000000.0f;


		sigma += resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha *
			(
			- sigmaEbackground_death / (2.f* param->Background_Sigma) + nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Background_Sigma )
			+     sigmaEcircle_death / (2.f* param->Shape_Sigma )     - nbNewPixels_death * 0.5f* log   ( 2 * CUDART_PI_F * param->Shape_Sigma      )
			);

	
	   if (isInside && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+nb_Cells]))
       {
			//add new one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);


						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)	//during the remove step, I already update the cell that are common to both ellipse
							{
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_new;
							}
						}
					}
				}
			}   		
			
			result_Ellipses->dev_axis0[index_obj] = newAxis0; 
			result_Ellipses->dev_axis1[index_obj] = newAxis1; 

	   }else{

		 //add old one back
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							 nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if ( nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)* nb_intersection] = area_old;
								dev_nbObjectsAtPixel_i[tmp1] =  nb_intersection+1;
							}
						}
					}
				}
			}
	   }
	}
}





__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void EnlargeShrink_on_Grid_bhat(const int* val, const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border,//input
							  Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 								//output
							   const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found, const unsigned int nb_Cells)											//const 
{
 
	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell
	const int index_obj = result_Ellipses->dev_index_not_available_spots[idx];

	if (idx <  nb_Cells && index_obj >= 0)
	{

		//--------------------------------get ellipse----------------------------------
		const float CenterX	= result_Ellipses->dev_x[index_obj];
		const float CenterY	= result_Ellipses->dev_y[index_obj];
		const float angle		= result_Ellipses->dev_angle[index_obj];
		const float axis0		= result_Ellipses->dev_axis0[index_obj]; 
		const float axis1		= result_Ellipses->dev_axis1[index_obj]; 
		//---------------------------------------------------------------------------------
		

		const unsigned int  global_y_origin =  (CenterY+border)/size_cell;	    				//index based on grid reference
		const unsigned int  global_x_origin =  (CenterX+border)/size_cell;	//index based on grid reference

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell

		if ( global_x_origin*size_cell >= w || global_y_origin*size_cell >= h) return;

		float newAxis0; float newAxis1;

		//1sd method (uniform ?)
		newAxis0 = param->Ellipse_MinWidth  + dev_rando[idx]				  * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth); 
		newAxis1 = param->Ellipse_MinHeight + dev_rando[idx+ nb_Cells] * ( param->Ellipse_MaxHeight - param->Ellipse_MinHeight);
		float sigma = newAxis0;//use as temporary
		newAxis0 = max(sigma, newAxis1); //axis 0 is bigger than axis 1
		newAxis1 = min(sigma, newAxis1);
		sigma = 0;

		//Cost to kill the current shape
		 unsigned int area_old = 0;unsigned int area_new = 0;
		 float resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0; 
		 
		//-----------------------Compute cost Interaction (resInter) and Data Cost (sigmaEbackground, sigmaEcircle)----------------------
		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;//must be signed
		const int border_boundary = param->Bhatt_Thickness;
		const int searchAreaExt = searchArea + border_boundary;

		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);

		float ratio0;float ratio1;float coordX; float coordY;
		int tmp1;
		int tmp2;unsigned int val_area;
		int current_val;
		unsigned int nb_intersection;
		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_old++;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area_new++;
			}
		}

		bool in_small, in_big, in_space, in_image;
		int position_area;
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{	  
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				//cost kill + remove
				coordX = CenterX + x;
				if ( coordX < 0 )  coordX = (int)(coordX - 1);
				else               coordX = (int)(coordX + 1);

				coordY = CenterY + y;
				if ( coordY < 0 )  coordY = (int)(coordY - 1);
				else               coordY = (int)(coordY + 1);

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
				in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
				in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
				in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

				if (in_small && in_space)//in ellipse and image with border
				{
					tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
					nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

					if ( in_image)// and in image
					{
						tmp1 = coordX+(w-2*border)*(coordY);
						current_val = val[tmp1];//origin at (0,0) image
						nbNewPixelsm1++;
						sigmaEcirclem1      += current_val * current_val;
						meanEcirclem1		+= current_val;
					}


					if ( nb_intersection > 1)//	(1 is itself)
					{
						 position_area = -1;
						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							if (val_area == area_old && position_area < 0)
							{
								position_area = index2;
							}else{
								resInter += (1.0f/ min ( area_old , val_area));
							}
						}
					}

					//remove at the same time
					if (nb_intersection > 0 )
					{
						if (nb_intersection > 1 && position_area >= 0)
							dev_minArea[tmp2+(w*h)*position_area] = dev_minArea[tmp2+(w*h)*( nb_intersection-1)];//swap last one with this one

						dev_nbObjectsAtPixel_i[tmp2] = nb_intersection-1;
					}

				}

				if (in_big && in_image && !in_small)//in ellipse and image with border
				{
					tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
					current_val = val[tmp1];
					nbNewPixelsm2++;
					sigmaEcirclem2      += current_val * current_val;
					meanEcirclem2		+= current_val;						
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma = -getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}
		//-------------------------------------------end cost to kill previous

		resInter = 0;
		nbNewPixelsm1 = 0;meanEcirclem1 = 0; sigmaEcirclem1 = 0;
		nbNewPixelsm2 = 0;meanEcirclem2 = 0; sigmaEcirclem2 = 0; 
		//-------------------------------------------cost to create new
		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				coordX = CenterX + x;
				if ( coordX < 0 )  coordX = (int)(coordX - 1);
				else               coordX = (int)(coordX + 1);

				coordY = CenterY + y;
				if ( coordY < 0 )  coordY = (int)(coordY - 1);
				else               coordY = (int)(coordY + 1);

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);
				in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f+border_boundary);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f+border_boundary);
				in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

				in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
				in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;

				if (in_small && in_space)//in ellipse and image with border
				{
					tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
					nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

					if ( in_image)// and in image
					{
						tmp1 = coordX+(w-2*border)*(coordY);
						current_val = val[tmp1];//origin at (0,0) image
						nbNewPixelsm1++;
						sigmaEcirclem1      += current_val * current_val;
						meanEcirclem1		+= current_val;
					}

					for (int index2 = 0;index2 < nb_intersection ; index2++)
					{	
						val_area = dev_minArea[tmp2+(w*h)*index2];
						resInter += (1.0f/ min ( area_new , val_area));
					} 

				}

				if (in_big && in_image && !in_small)//in ellipse and image with border
				{
					tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
					current_val = val[tmp1];
					nbNewPixelsm2++;
					sigmaEcirclem2      += current_val * current_val;
					meanEcirclem2		+= current_val;						
				}
			}
		}

		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			sigma += getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);			
		}

	
	   if (nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma, temperature , 1 , dev_rando[idx+2*nb_Cells]))
       {
			  
			//add new one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (newAxis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (newAxis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);



						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_new;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}   
	
			result_Ellipses->dev_axis0[index_obj] = newAxis0; 
			result_Ellipses->dev_axis1[index_obj] = newAxis1; 
		}
		

		else{
			
	 		//add old one
			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = CenterX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = CenterY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border and in the new one only (!in_both)
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area_old;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}
						
						
					}
				}
			}
		}		
	}
}


//the new one (not used)
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Birth_on_Grid_RJMCMC_bhat_nonHo(const int* val,const int w,const int h,const float *dev_rando, const int size_cell,const int offsetw,const int offseth,const int border, const unsigned int  nb_Cells, //input
				   Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 																							//output
				   const unsigned int* dev_array_quadtree_x,const unsigned int* dev_array_quadtree_y,const unsigned int* dev_array_quadtree_w, const unsigned int* dev_array_quadtree_h, const unsigned int nb_thread_total, const unsigned int offset,
				    const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found)																								//const 
{

	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x) ;//idx is the index of the first pixel of the grill cell


	if (idx < nb_thread_total)
	{
		const unsigned int  global_y_origin  = dev_array_quadtree_y[idx + offset] ;			//index based on space reference
		const unsigned int  global_x_origin  = dev_array_quadtree_x[idx + offset] ;			//index based on space reference

		const unsigned int local_size_cell_w = dev_array_quadtree_w[idx + offset];						//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = dev_array_quadtree_h[idx + offset];						//real sizeh in pixel of the current cell


		float centerX = dev_rando[idx];
		float centerY = dev_rando[idx + nb_thread_total];

		int global_y_position =  local_size_cell_h * centerY + global_y_origin - border;		
		int global_x_position =  local_size_cell_w * centerX + global_x_origin - border;	   

		//--------------------------------create ellipse----------------------------------
		centerX += global_x_position;
		centerY += global_y_position;
		float angle   = dev_rando[idx+ (2*nb_thread_total)]* 2 * CUDART_PI_F ;
		float axis0   = param->Ellipse_MinWidth  + dev_rando[idx+ (3*nb_thread_total)] * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth);
		float axis1   = param->Ellipse_MinHeight + dev_rando[idx+ (4*nb_thread_total)] * ( param->Ellipse_MaxHeight - param->Ellipse_MinHeight);

		float tmp_f;
		if (axis0 < axis1)//swap, axis0 must be bigger than axis1
		{
			tmp_f = axis0;axis0 = axis1;axis1 = tmp_f;
		}


		float resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0;

		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;
		const int border_boundary = param->Bhatt_Thickness;

		const int searchAreaExt = searchArea + border_boundary;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;

		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int tmp1;
		int tmp2;
		unsigned int val_area;
		int current_val;//need to be signed because we substract it to something else

		unsigned int nb_intersection;//int other_area;
		bool in_small, in_big, in_space, in_image;

		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt + 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				{
					coordX = centerX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = centerY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;


					if (in_small && in_space)//in ellipse and image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// in image
						{
							tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
							current_val = val[tmp1];
							nbNewPixelsm1++;
							sigmaEcirclem1      += current_val * current_val;
							meanEcirclem1		+= current_val;
						}

						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += __frcp_rn((int)min ( area , val_area));//(1.0f/ min ( area , val_area));
						}							
					}

					if (in_big && in_image&& !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		float sigma_birth = param->Intersections_Alpha + param->Bhatt_Alpha;
		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			
			sigma_birth = getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);

		}

		unsigned int indexGlobalGrid = int(centerX)/size_cell + ceil(float(w)/size_cell)*(int(centerY)/size_cell);


		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[indexGlobalGrid];

		if ( nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma_birth , temperature , float(local_size_cell_w*local_size_cell_h)/(nb_Ellipse_in_this_cell+1), dev_rando[idx+5*nb_thread_total]))
		{
			unsigned int index_spot_available = result_Ellipses->dev_index_available_spots[idx];

			//save the new ellipse
			result_Ellipses->dev_x[index_spot_available]				= centerX;	 
			result_Ellipses->dev_y[index_spot_available]				= centerY;
			result_Ellipses->dev_angle[index_spot_available]			= angle;
			result_Ellipses->dev_axis0[index_spot_available]			= axis0;
			result_Ellipses->dev_axis1[index_spot_available]			= axis1;

			result_Ellipses->dev_available_spot[index_spot_available]	= false;

			dev_nbObjectsAtCell_i[indexGlobalGrid] = nb_Ellipse_in_this_cell + 1;


			unsigned int nb_overlap;

			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = centerX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = centerY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_overlap = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_overlap+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_overlap] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_overlap+1;
							}
						}

					}
				}	
			}   
		}
	}
}



	 

//this one is "optimize" to go faster by reading less global memory
//not cut
//size_cell is the smallest cell possible (used as reference), whereas size_cell_qt is the current real size depending of witch level from the qt we work on.
__global__ //__launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Birth_on_Grid_RJMCMC_bhat_nonHo(const int* val,const int w,const int h,const float *dev_rando, const int size_cell,const int offsetw,const int offseth,const int border, const unsigned int  nb_Cells, //input
				   Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 																							//output
				   const unsigned int* index_grid,const unsigned int size_cell_qt, const unsigned int nb_thread_total, const unsigned int offset, const float ratio_birth_death,
				    const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found)																							//const 
{

	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x);//idx is the index of the first pixel of the grill cell

	if (idx < nb_thread_total )
	{
		const unsigned int pos_pixel = index_grid[idx + offset];

		//!!!! w and h are already with 2*border !!!!
		const unsigned int  global_y_origin  = pos_pixel / w ;								
		const unsigned int  global_x_origin  = pos_pixel - global_y_origin*w ;	
		
		const unsigned int local_size_cell_w = min(size_cell_qt,abs(int(w-global_x_origin) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell_qt,abs(int(h-global_y_origin) ) );		//real sizeh in pixel of the current cell


		const int global_y_position =  local_size_cell_h * dev_rando[idx];																//a pixel in this cell, local origin (x)
		const int global_x_position =  local_size_cell_w * dev_rando[idx + nb_thread_total];	

		//--------------------------------create Ellipse----------------------------------
		const float centerX = global_x_origin + global_x_position - border;//in image coordinate
		const float centerY = global_y_origin + global_y_position - border;//in image coordinate

	if ( centerX >= w-2*border || centerY >= h-2*border || centerX < -border || centerY < -border) return;

		float angle;
		if (ELLIPSE_STEP_ANGLE > 0)
			angle   = int(dev_rando[idx+ 2*nb_thread_total] * ELLIPSE_STEP_ANGLE) * CUDART_PI_F / ELLIPSE_STEP_ANGLE ;
		else
			angle   = dev_rando[idx+ 2*nb_thread_total] * CUDART_PI_F ;

		float axis0   = param->Ellipse_MinWidth  + dev_rando[idx+ 3*nb_thread_total] * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth);
		float axis1   = param->Ellipse_MinHeight + dev_rando[idx+ 4*nb_thread_total] * ( param->Ellipse_MaxHeight - param->Ellipse_MinHeight);

		float resInter;
		if (axis0 < axis1)//swap, axis0 must be bigger than axis1
		{
			resInter = axis0;axis0 = axis1;axis1 = resInter;
		}


		resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0;

		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;
		const int border_boundary = param->Bhatt_Thickness;

		const int searchAreaExt = searchArea + border_boundary;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);

		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;

		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int tmp1;
		int tmp2;
		unsigned int val_area;
		int current_val;//need to be signed because we substract it to something else

		unsigned int nb_intersection;//int other_area;
		bool in_small, in_big, in_space, in_image;

		for ( int y = -searchAreaExt - 0.5f ; y <= searchAreaExt + 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt + 0.5f; x++ )
			{
				{
					coordX = centerX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = centerY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f + border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f + border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;


					if (in_small && in_space)//in ellipse and image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// in image
						{
							tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
							current_val = val[tmp1];
							nbNewPixelsm1++;
							sigmaEcirclem1      += current_val * current_val;
							meanEcirclem1		+= current_val;
						}

						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area , val_area));
						}	
					}

					if (in_big && in_image&& !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		float sigma_birth = param->Intersections_Alpha + param->Bhatt_Alpha;
		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			
			sigma_birth = getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);
		}

		if ( resInter > param->Hard_Core ) return;

		tmp1 = global_x_origin + ceil(float(w)/size_cell)*global_y_origin;
		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[tmp1];
		

		if ( nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma_birth , temperature , float(local_size_cell_w*local_size_cell_h)/(nb_Ellipse_in_this_cell+1), dev_rando[idx+5*nb_thread_total]))
		{
			unsigned int index_spot_available = result_Ellipses->dev_index_available_spots[idx];

			//save the new ellipse
			result_Ellipses->dev_x[index_spot_available]				= centerX;	 
			result_Ellipses->dev_y[index_spot_available]				= centerY;
			result_Ellipses->dev_angle[index_spot_available]			= angle;
			result_Ellipses->dev_axis0[index_spot_available]			= axis0;
			result_Ellipses->dev_axis1[index_spot_available]			= axis1;
			result_Ellipses->dev_available_spot[index_spot_available]	= false;

			dev_nbObjectsAtCell_i[tmp1] = nb_Ellipse_in_this_cell + 1;

			for ( int y = -searchArea - 0.5f; y <= searchArea + 0.5f; y++ )
			{
				for ( int x = -searchArea - 0.5f; x <= searchArea+ 0.5f; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = centerX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = centerY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}

					}
				}	
			}   
		}
	}
}



//the new one
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Birth_on_Grid_RJMCMC_bhat(const int* val,const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border, const unsigned int  nb_Cells, //input
				   Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i,																						//output
				   const float temperature, bhatt_Parameter *param, const unsigned int nb_Objects_found)																								//const 
{

	const unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;//idx is the index of the first pixel of the grill cell

	if (idx <= nb_Cells )
	{
		//get location
		unsigned int  global_y_origin =  (unsigned int(__fdiv_rn(idx , nb_cell_w)) << 1);	    				//index based on grid reference
		unsigned int  global_x_origin = (idx << 1) - nb_cell_w*global_y_origin + offsetw;	//index based on grid reference
		global_y_origin+= offseth;

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell


		float centerX = dev_rando[idx];
		float centerY = dev_rando[idx + nb_Cells];

		int global_y_position =  local_size_cell_h * centerY;																//a pixel in this cell, local origin (x)
		int global_x_position =  local_size_cell_w * centerX;																//a pixel in this cell, local origin (y)
		global_x_position+=global_x_origin*size_cell - border;global_y_position+=global_y_origin*size_cell - border;		//a pixel in this cell, global (image) origin

		//--------------------------------create ellipse----------------------------------
		centerX += global_x_position;//in image coordinate
		centerY += global_y_position;//in image coordinate


		if ( centerX >= w-2*border || centerY >= h-2*border || centerX < -border || centerY < -border) return;

		float angle;
		if (ELLIPSE_STEP_ANGLE > 0)
			angle   = int(dev_rando[idx+ 2*nb_Cells] * ELLIPSE_STEP_ANGLE) * CUDART_PI_F / ELLIPSE_STEP_ANGLE ;
		else
			angle   = dev_rando[idx+ 2*nb_Cells] * CUDART_PI_F ;

		float axis0   = param->Ellipse_MinWidth  + dev_rando[idx+ 3*nb_Cells] * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth);
		
		float axis1 = param->Ellipse_MaxHeight;

		axis1   = param->Ellipse_MinHeight + dev_rando[idx+ 4*nb_Cells] * ( axis1 - param->Ellipse_MinHeight);

		float tmp_f;
		if (axis0 < axis1)//swap, axis0 must be bigger than axis1
		{
			tmp_f = axis0;axis0 = axis1;axis1 = tmp_f;
		}


		float resInter = 0;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0;

		const int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;
		const int border_boundary = param->Bhatt_Thickness;

		const int searchAreaExt = searchArea + border_boundary;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);
		//float axis1vec[2];axis0vec[1]=axis0vec[1];axis1vec[1]=-axis0vec[0];
		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;

		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int tmp1;
		int tmp2;
		unsigned int val_area;
		int current_val;//need to be signed because we substract it to something else

		unsigned int nb_intersection;//int other_area;
		bool in_small, in_big, in_space, in_image;

		for ( int y = -searchAreaExt - 0.5f; y <= searchAreaExt+ 0.5f; y++ )
		{
			for ( int x = -searchAreaExt - 0.5f; x <= searchAreaExt+ 0.5f; x++ )
			{
				{
					coordX = centerX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = centerY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);
					in_small = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f+border_boundary);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f+border_boundary);
					in_big = ratio0 * ratio0 + ratio1 * ratio1 <= 1.f;

					in_image = coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border;
					in_space = coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border;


					if (in_small && in_space)//in ellipse and image with border
					{
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( in_image)// in image
						{
							tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
							current_val = val[tmp1];
							nbNewPixelsm1++;
							sigmaEcirclem1      += current_val * current_val;
							meanEcirclem1		+= current_val;
						}

						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += (1.0f/ min ( area , val_area));
						}	
					}

					if (in_big && in_image&& !in_small)//in ellipse and image with border
					{
						tmp1 = coordX+(w-2*border)*coordY;//origin at (0,0) image
						current_val = val[tmp1];
						nbNewPixelsm2++;
						sigmaEcirclem2      += current_val * current_val;
						meanEcirclem2		+= current_val;						
					}
				}
			}
		}

		float sigma_birth = param->Intersections_Alpha + param->Bhatt_Alpha;
		if	(nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS) 
		{
			meanEcirclem2/=nbNewPixelsm2;	meanEcirclem1/=nbNewPixelsm1;
			sigmaEcirclem1 = (float(nbNewPixelsm1-1)/nbNewPixelsm1)*(sigmaEcirclem1/nbNewPixelsm1 - meanEcirclem1*meanEcirclem1);
			sigmaEcirclem2 = (float(nbNewPixelsm2-1)/nbNewPixelsm2)*(sigmaEcirclem2/nbNewPixelsm2 - meanEcirclem2*meanEcirclem2);
			if ( sigmaEcirclem1 < 0.0000001f) sigmaEcirclem1  = 0.0000001f;
			if ( sigmaEcirclem2 < 0.0000001f) sigmaEcirclem2  = 0.0000001f;

			if ( resInter > param->Hard_Core )resInter = 1.0f;

			
			sigma_birth = getMesure( meanEcirclem1, meanEcirclem2, sigmaEcirclem1, sigmaEcirclem2, resInter, param);
		}

		if ( resInter > param->Hard_Core ) return;

		tmp1 = global_x_origin + 2*nb_cell_w*global_y_origin;
		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[tmp1];

		if ( nbNewPixelsm1 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && evaluateGreen( sigma_birth , temperature , float(local_size_cell_w*local_size_cell_h)/(nb_Ellipse_in_this_cell+1), dev_rando[idx+5*nb_Cells]))
		{
			unsigned int index_spot_available = result_Ellipses->dev_index_available_spots[idx];

			//save the new ellipse
			result_Ellipses->dev_x[index_spot_available]				= centerX;	 
			result_Ellipses->dev_y[index_spot_available]				= centerY;
			result_Ellipses->dev_angle[index_spot_available]			= angle;
			result_Ellipses->dev_axis0[index_spot_available]			= axis0;
			result_Ellipses->dev_axis1[index_spot_available]			= axis1;

			result_Ellipses->dev_available_spot[index_spot_available]	= false;

			dev_nbObjectsAtCell_i[tmp1] = nb_Ellipse_in_this_cell + 1;

			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = centerX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = centerY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_intersection = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_intersection+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_intersection] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_intersection+1;
							}
						}

					}
				}	
			}   
		}
	}
}




//the new one
__global__// __launch_bounds__(BOUND_MAX_THREADS,BOUND_MIN_BLOCK_PER_SM)
void Birth_on_Grid_RJMCMC_bayesian(const int* val,const int w,const int h,const float *dev_rando,const int nb_cell_w,const int nb_cell_h,const int size_cell,const int offsetw,const int offseth,const int border, const unsigned int  nb_Cells, //input
				   Ellipse_result *result_Ellipses, unsigned int* dev_minArea, unsigned int* dev_nbObjectsAtPixel_i,unsigned int* dev_nbObjectsAtCell_i, 																							//output
				    const float temperature, bayesian_Parameter *param, const unsigned int nb_Objects_found)																								//const 
{
	const unsigned int idx = threadIdx.x + (blockIdx.x*blockDim.x) ;//idx is the index of the first pixel of the grill cell


	if (idx < nb_Cells)
	{
		//get location
		unsigned int  global_y_origin =  (unsigned int(__fdiv_rn(idx , nb_cell_w)) << 1);	    				//index based on grid reference
		unsigned int  global_x_origin = (idx << 1) - nb_cell_w*global_y_origin + offsetw;	//index based on grid reference
		global_y_origin+= offseth;

		const unsigned int local_size_cell_w = min(size_cell,abs(w-int(global_x_origin*size_cell) ) );		//real sizew in pixel of the current cell
		const unsigned int local_size_cell_h = min(size_cell,abs(h-int(global_y_origin*size_cell) ) );		//real sizeh in pixel of the current cell


		float centerX = dev_rando[idx];
		float centerY = dev_rando[idx + nb_Cells];

		int global_y_position =  local_size_cell_h * centerY;																//a pixel in this cell, local origin (x)
		int global_x_position =  local_size_cell_w * centerX;																//a pixel in this cell, local origin (y)
		global_x_position+=global_x_origin*size_cell - border;global_y_position+=global_y_origin*size_cell - border;		//a pixel in this cell, global (image) origin

		//--------------------------------create ellipse----------------------------------
		centerX += global_x_position;//in image coordinate
		centerY += global_y_position;//in image coordinate
		float angle   = dev_rando[idx+ (2*nb_Cells)]* CUDART_PI_F ;
		float axis0   = param->Ellipse_MinWidth  + dev_rando[idx+ (3*nb_Cells)] * ( param->Ellipse_MaxWidth  - param->Ellipse_MinWidth);
		float axis1   = param->Ellipse_MinHeight + dev_rando[idx+ (4*nb_Cells)] * ( param->Ellipse_MaxHeight - param->Ellipse_MinHeight);

		float tmp_f;
		if (axis0 < axis1)//swap, axis0 must be bigger than axis1
		{
			tmp_f = axis0;axis0 = axis1;axis1 = tmp_f;
		}


		float resInter = 0;
		unsigned int nbNewPixels_birth = 0;float sigmaEbackground_birth = 0; float sigmaEcircle_birth = 0;

		int searchArea = max(param->Ellipse_MaxWidth, param->Ellipse_MaxHeight) + 0.5f;
		float axis0vec[2];axis0vec[0]=cosf(angle);axis0vec[1]=sinf(angle);

		float ratio0;float ratio1;float coordX; float coordY;
		unsigned int area = 0;

		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
				ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

				if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f) area++;
			}
		}

		int tmp1;
		int tmp2;
		unsigned int val_area;
		int current_val;//need to be signed because we substract it to something else

		unsigned int nb_intersection;//int other_area;

		for ( int y = -searchArea ; y <= searchArea; y++ )
		{
			for ( int x = -searchArea ; x <= searchArea; x++ )
			{
				{
					coordX = centerX + x;
					if ( coordX < 0 )  coordX = (int)(coordX - 1);
					else               coordX = (int)(coordX + 1);

					coordY = centerY + y;
					if ( coordY < 0 )  coordY = (int)(coordY - 1);
					else               coordY = (int)(coordY + 1);

					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);


					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1.f&& coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
					{

						tmp1 = coordX+(w-(int)(2*border))*(coordY);//origin at (0,0) image
						tmp2 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
						nb_intersection = dev_nbObjectsAtPixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < w - 2*border && coordY < h - 2*border)// in image
						{
							if (nb_intersection == 0)
							{
								current_val = val[tmp1];
								nbNewPixels_birth++;
								sigmaEbackground_birth  += (current_val - param->Background_Avg)	*	(current_val - param->Background_Avg) ;
								sigmaEcircle_birth      += (current_val - param->Shape_Avg)			*	(current_val - param->Shape_Avg)     ;
							}
						}


						for (int index2 = 0;index2 < nb_intersection ; index2++)
						{	
							val_area = dev_minArea[tmp2+(w*h)*index2];
							resInter += __frcp_rn((int)min ( area , val_area));//(1.0f/ min ( area , val_area));
						}							


					}
				}
			}
		}

		if ( resInter > param->Hard_Core )resInter = 10000000000.0f;


		float sigma_birth = resInter*param->Intersections_Alpha 
			+ param->Bayesian_Alpha *
			(
			- sigmaEbackground_birth / (2.f* param->Background_Sigma) + nbNewPixels_birth * 0.5f* log   ( 2 * CUDART_PI_F * param->Background_Sigma )
			+     sigmaEcircle_birth / (2.f* param->Shape_Sigma )     - nbNewPixels_birth * 0.5f* log   ( 2 * CUDART_PI_F * param->Shape_Sigma      )
			);

		tmp1 = 	global_x_origin + 2*nb_cell_w*global_y_origin;
		unsigned int nb_Ellipse_in_this_cell = dev_nbObjectsAtCell_i[tmp1];


		unsigned int index_spot_available = result_Ellipses->dev_index_available_spots[idx];
		if ( evaluateGreen( sigma_birth , temperature , float(local_size_cell_w*local_size_cell_h)/(nb_Ellipse_in_this_cell+1), dev_rando[idx+5*nb_Cells]))
		{
			//save the new ellipse
			result_Ellipses->dev_x[index_spot_available]				= centerX;	 
			result_Ellipses->dev_y[index_spot_available]				= centerY;
			result_Ellipses->dev_angle[index_spot_available]			= angle;
			result_Ellipses->dev_axis0[index_spot_available]			= axis0;
			result_Ellipses->dev_axis1[index_spot_available]			= axis1;

			result_Ellipses->dev_available_spot[index_spot_available]	= false;

			
			dev_nbObjectsAtCell_i[tmp1] = nb_Ellipse_in_this_cell + 1;

			unsigned int nb_overlap;

			for ( int y = -searchArea ; y <= searchArea; y++ )
			{
				for ( int x = -searchArea ; x <= searchArea; x++ )
				{
					ratio0 = (axis0vec[0]*x + axis0vec[1]*y) / (axis0 + 0.5f);
					ratio1 = (axis0vec[1]*x - axis0vec[0]*y) / (axis1 + 0.5f);

					if ( ratio0 * ratio0 + ratio1 * ratio1 <= 1.f)//that is in_ellipse
					{
						coordX = centerX + x;
						if ( coordX < 0 )  coordX = (int)(coordX - 1);
						else               coordX = (int)(coordX + 1);

						coordY = centerY + y;
						if ( coordY < 0 )  coordY = (int)(coordY - 1);
						else               coordY = (int)(coordY + 1);

						if (coordX >= -border && coordX < w - border && coordY >= -border && coordY < h - border)//in ellipse and image with border
						{
							tmp1 = coordX+border+w*(coordY+border);//origin at (0,0) space (with border) !
							nb_overlap = dev_nbObjectsAtPixel_i[tmp1];

							if (nb_overlap+1 < NB_OVERLAP_MAX_ON_PIXEL)
							{
								dev_minArea[tmp1+(w*h)*nb_overlap] = area;
								dev_nbObjectsAtPixel_i[tmp1] = nb_overlap+1;
							}
						}

					}
				}	
			}   
		}
	}
}						  