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


//-----------------------------------------------------------------------------------------------//
//-------------------------------------GLOBAL DECLARATION----------------------------------------//
//-----------------------------------------------------------------------------------------------//
//------------On Host  
static bool isInit = false;
static METHOD_DATATERM method_dterm;
curandGenerator_t gen;
static cudaDeviceProp host_prop;
static curandStatus status;
static cudaEvent_t start_event;
static cudaEvent_t stop_event;
static unsigned int host_capacity;
unsigned int* host_index_spots_available;
bool* host_spots_available;
Const_Parameter* host_param;
static quadtree *host_quadtree = NULL;	//for image or tree3D
float ratio_birth_death;

unsigned int nb_Objects_found;
unsigned int max_size_global_result;
int nb_Cells;int size_cell;
static const int grid[][2] =  {{0,0},{0,1},{1,0},{1,1}}; int nb_cell_grid[4][3];
int nb_grid_w;int nb_grid_h;
int host_width; int host_height; int host_depth;int host_border;
int host_nb_thread_per_block;//depend on the arch

//for image (ellipse)
int* host_image;

int host_sizeArray;
int host_origin_x, host_origin_y,  host_origin_z;


Ellipse_result *host_result_Ellipses = NULL;//global result
//---------------------------------

//-----------On Device
float* dev_rand;
float* dev_BirthMap;
unsigned int* dev_minArea;
float* dev_dataAtPixel;
unsigned int* dev_nbObjectsAtPixel_i;
unsigned int* dev_nbObjectsAtCell_i;
Const_Parameter *dev_param;

//for images
int* dev_image;

float *dev_arrayX, *dev_arrayY, *dev_arrayZ;
unsigned int *dev_array_quadtree_x, *dev_array_quadtree_y, *dev_array_quadtree_w, *dev_array_quadtree_h;
float* manifoldTexture = NULL;


Ellipse_result *dev_result_Ellipses = NULL;//global result
//-----------------------------------



//-----------------------------------------------------------------------------------------------//
//---------------------------------------END DECLARATION-----------------------------------------//
//-----------------------------------------------------------------------------------------------//


void checkCUDAError(const char *msg)
{
	cudaThreadSynchronize();
	cudaError_t err = cudaGetLastError();
	if( cudaSuccess != err) 
	{
		fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
		exit(EXIT_FAILURE);
	}                         
}


//return memory free, and set the number of thread_per_block
unsigned int initDevice(int iddevice)
{
	cudaDeviceReset ();
	cuInit(0); // Initialize CUDA
	//---------------------just for display purpose-----------------------------
	//int iddevice = 1;//1 quadro, 0 GTX 480			//<----------- select GPU here
	int max; cudaGetDeviceCount(&max);
	iddevice = iddevice > max - 1 ? max - 1 : iddevice ;
	printf("use device %d\n",iddevice);
	cudaSetDevice(iddevice);
	CUdevice device;
	CUcontext context;
	cuDeviceGet( &device, iddevice ); // Get handle for device
	cuCtxCreate( &context, 0, device ); // Create context

	unsigned int free2 = 0;
	unsigned int total;
	cuMemGetInfo((size_t*)&free2,(size_t*) &total);

    int deviceIndex;
    cudaGetDevice(&deviceIndex);

    cudaGetDeviceProperties(&host_prop, deviceIndex);
	host_capacity = (host_prop.major*10+host_prop.minor)*10;
	printf("capacity:\t%.1f\n",host_capacity/10.);
	printf("number of SM:\t%d\n",host_prop.multiProcessorCount);

	//like in MACRO_Definition for __bound__ (see GPU_MACRO_Definition.h)
	if (host_capacity >=200)
		host_nb_thread_per_block = NB_THEARD_PER_BLOCK ;
	else if	(host_capacity >=120)
		host_nb_thread_per_block = NB_THEARD_PER_BLOCK /2;
	else
		host_nb_thread_per_block = NB_THEARD_PER_BLOCK /4;
	//---------------------------------------------------------------------------

	return free2;
}

void prepareEllipsesCandidateForGrid_i(int grid2run, int w, int h)
{
	checkCUDAError("error at prepare");
	unsigned int* nbResultperCell = (unsigned int*) malloc(4*nb_Cells*sizeof(unsigned int));
	int* objects2pick = (int*) malloc(4*nb_Cells*sizeof(int));

	int* index_spots_not_available = (int*) malloc(nb_Cells*sizeof(int));
	memset(index_spots_not_available, -1, nb_Cells*sizeof(int));

	cudaMemcpy(host_spots_available, host_result_Ellipses->dev_available_spot,max_size_global_result*sizeof(bool),cudaMemcpyDeviceToHost);
	cudaMemcpy(nbResultperCell,	dev_nbObjectsAtCell_i ,4*nb_Cells*sizeof(unsigned int),cudaMemcpyDeviceToHost);

	for (int i=0;i<4*nb_Cells;i++)
		objects2pick[i] = rand()%(1+nbResultperCell[i]);//pick a number between [0, nb_object_in_cell_i], this is the selected object

	unsigned int index_cell_object;
	unsigned int index_grid;

	float* centerX = (float*) malloc(max_size_global_result*sizeof(float));
	float* centerY = (float*) malloc(max_size_global_result*sizeof(float));
	cudaMemcpy(centerX, host_result_Ellipses->dev_x,max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);
	cudaMemcpy(centerY, host_result_Ellipses->dev_y,max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);

	int idx, idy;
	int count = 0;
	for (unsigned int i=0;i<max_size_global_result;i++)
	{
		if (!host_spots_available[i])//used
		{

			idx = (centerX[i]+host_border)/size_cell;
			idy = (centerY[i]+host_border)/size_cell;

			index_grid = 2*(idx&1) + idy&1;

			index_cell_object =	idx + idy*2*grid[grid2run][0];
			if(index_grid == grid2run)//same grid [0,3]
			{
				if (objects2pick[index_cell_object] == 0)
					index_spots_not_available[count++] = i;

				objects2pick[index_cell_object]--;
			}
		}
	}

	free(centerX);free(centerY);
	cudaMemcpy(host_result_Ellipses->dev_index_not_available_spots, index_spots_not_available, nb_Cells*sizeof(int), cudaMemcpyHostToDevice);//-1 means no ellipse
	free(objects2pick);

	free(index_spots_not_available);
	free(nbResultperCell);//free(spots_available);


}


float getEnergy(float centerX, float centerY, float angle,float axis0,float axis1, unsigned int* Pixel_i, unsigned int* minArea, bool verbose, bool onlyDataTerm = false)
{
	//host_image
	int searchArea = 0;
	
	if (method_dterm == METHOD_DATATERM::BAYESIAN) 
		searchArea = max(((bayesian_Parameter*)host_param)->Ellipse_MaxWidth, ((bayesian_Parameter*)host_param)->Ellipse_MaxHeight) + 0.5;//keep it signed it, if not the for loop (int y, int x) above will failed
	
	if (method_dterm == METHOD_DATATERM::BHATTACHARYA) 
		searchArea = max( ((bhatt_Parameter*)host_param)->Ellipse_MaxWidth,  ((bhatt_Parameter*)host_param)->Ellipse_MaxHeight) + 0.5;//keep it signed it, if not the for loop (int y, int x) above will failed
	
	float axis0vec[2];axis0vec[0]=cos(angle);axis0vec[1]=sin (angle);
	float axis1vec[2];axis1vec[0]=axis0vec[1];axis1vec[1]=-axis0vec[0];
	float ratio0;float ratio1;float coordX; float coordY;
	int tmp1;int tmp2;	unsigned int  nb_intersection;
	float sigmaEbackground_birth = 0;
	float sigmaEcircle_birth = 0;
	unsigned int nbNewPixels_birth = 0;
	int current_val;
	unsigned int area = 0;
	float resInter = 0;
	unsigned int val_area;

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

				ratio0 /=   (axis0 + 0.5f );
				ratio1 /=   (axis1 + 0.5f );

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


	if (method_dterm == METHOD_DATATERM::BAYESIAN) 
	{

		for ( int y = -searchArea - 0.5f; y <= searchArea + 0.5f; y++ )
		{
			for ( int x = -searchArea - 0.5f; x <= searchArea + 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.5 );
					ratio1 = (axis1vec[0]*x + axis1vec[1]*y) / (axis1 + 0.5 );

					if (ratio0 * ratio0 + ratio1 * ratio1 <= 1. && coordX >= -host_border && coordX < host_width + host_border && coordY >= -host_border && coordY < host_height + host_border)//in ellipse and image with border
					{
						tmp1 = coordX+host_width*coordY;//origin at (0,0) image
						tmp2 = coordX+host_border+(host_width+2*host_border)*(coordY+host_border);//origin at (0,0) space (with border) !
						nb_intersection = Pixel_i[tmp2];

						if ( coordX >= 0 && coordY >= 0 && coordX < host_width && coordY < host_height)// in image
						{
							if (nb_intersection == 1)   //only itself
							{
								current_val = host_image[tmp1];
								nbNewPixels_birth++;
								sigmaEbackground_birth  += (current_val - ((bayesian_Parameter*)host_param)->Background_Avg)	    *	(current_val - ((bayesian_Parameter*)host_param)->Background_Avg) ;
								sigmaEcircle_birth      += (current_val - ((bayesian_Parameter*)host_param)->Shape_Avg)			*	(current_val - ((bayesian_Parameter*)host_param)->Shape_Avg)     ;
							}
						}

						if ( nb_intersection > 1)//	(1 is itself)
						{
							bool already_jump = false;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = minArea[tmp2+(host_width+2*host_border)*(host_height+2*host_border)*index2];
								if (val_area == area && !already_jump)
								{
									already_jump = true;
								}else{
									resInter += (1.0 / min ( area , val_area));
								}
							}
						}				

					}
				}
			}
		}

		if ( resInter > host_param->Hard_Core )resInter = 10000000000.0;


		return resInter*host_param->Intersections_Alpha 
			+ ((bayesian_Parameter*)host_param)->Bayesian_Alpha *
			(
			- sigmaEbackground_birth / (2. * ((bayesian_Parameter*)host_param)->Background_Sigma) + nbNewPixels_birth * 0.5 * log   ( 2 * CUDART_PI_F * ((bayesian_Parameter*)host_param)->Background_Sigma )
			+     sigmaEcircle_birth / (2. * ((bayesian_Parameter*)host_param)->Shape_Sigma )     - nbNewPixels_birth * 0.5 * log   ( 2 * CUDART_PI_F * ((bayesian_Parameter*)host_param)->Shape_Sigma      )
			);
	}


	 float sigma_birth = 0;

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA) 
	{
		const int searchAreaExt = searchArea + ((bhatt_Parameter*)host_param)->Bhatt_Thickness;

		bool in_small, in_big, in_space, in_image;
		unsigned int nbNewPixelsm1 = 0;float meanEcirclem1 = 0; float sigmaEcirclem1 = 0;
		unsigned int nbNewPixelsm2 = 0;float meanEcirclem2 = 0; float sigmaEcirclem2 = 0;

		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.;

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

					in_image = coordX >= 0 && coordY >= 0 && coordX < host_width  && coordY < host_height;
					in_space = coordX >= -host_border && coordX < host_width + host_border && coordY >= -host_border && coordY < host_height + host_border;

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

						if ( in_image)// in image
						{
							tmp1 = coordX+host_width*coordY;//origin at (0,0) image
							current_val = host_image[tmp1];
							nbNewPixelsm1++;
							sigmaEcirclem1      += current_val * current_val;
							meanEcirclem1		+= current_val;
						}

						if ( nb_intersection > 1 && nb_intersection < NB_OVERLAP_MAX_ON_PIXEL)//	(1 is itself)
						{
							bool already_jump = false;
							//tmp1 = area;
							for (int index2 = 0;index2 < nb_intersection ; index2++)
							{	
								val_area = minArea[tmp2+(host_width+2*host_border)*(host_height+2*host_border)*index2];
								if (val_area == area && !already_jump)
								{
									already_jump = true;
								}else{
									resInter += (1.f / min ( area , val_area));
								}
							}
						}	
					}

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

		sigma_birth = ((bhatt_Parameter*)host_param)->Intersections_Alpha + ((bhatt_Parameter*)host_param)->Bhatt_Alpha;
		if	(nbNewPixelsm2 > ELLIPSE_MIN_NB_PIXELS && nbNewPixelsm1 > 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.0000001 ) sigmaEcirclem1  = 0.0000001;
			if ( sigmaEcirclem2 < 0.0000001 ) sigmaEcirclem2  = 0.0000001;

			if ( resInter > host_param->Hard_Core )resInter = 1.0;


	sigma_birth = 100.f*((meanEcirclem1-meanEcirclem2)*(meanEcirclem1-meanEcirclem2)/(4.f*(sigmaEcirclem1+sigmaEcirclem2))
		-0.5f*logf((2.f*sqrtf(sigmaEcirclem1*sigmaEcirclem2))/(sigmaEcirclem1+sigmaEcirclem2)));
	float d0 = ((bhatt_Parameter*)host_param)->Bhatt_Threshold;


	if ( 
		ELLIPSE_WRONG_DATA_CONTRAST(meanEcirclem1,meanEcirclem2) || ELLIPSE_WRONG_DATA_BOUND(meanEcirclem1,meanEcirclem2,sigmaEcirclem1,sigmaEcirclem2)
		) 
		if (onlyDataTerm)
			return ((bhatt_Parameter*)host_param)->Bhatt_Alpha;
		else
			return ((bhatt_Parameter*)host_param)->Intersections_Alpha + ((bhatt_Parameter*)host_param)->Bhatt_Alpha;
	else
		if (sigma_birth <= d0)
		{
			if (onlyDataTerm)
				return  ((bhatt_Parameter*)host_param)->Bhatt_Alpha*(1.f-sigma_birth / d0 ) ;
			else
				return  ((bhatt_Parameter*)host_param)->Bhatt_Alpha*(1.f-sigma_birth / d0 ) + ((bhatt_Parameter*)host_param)->Intersections_Alpha*resInter;			
		}
		else 
		{
			if (onlyDataTerm)
				return ((bhatt_Parameter*)host_param)->Bhatt_Alpha*(exp((d0-sigma_birth)/(d0)) - 1.f) ;
			else	
				return ((bhatt_Parameter*)host_param)->Bhatt_Alpha*(exp((d0-sigma_birth)/(d0)) - 1.f) + ((bhatt_Parameter*)host_param)->Intersections_Alpha*resInter;	 
		}


		}
		return sigma_birth;
	}


	return sigma_birth;
}

extern "C" unsigned int getNbResult()
{

	if	(host_result_Ellipses != NULL)
		cudaMemcpy(host_spots_available, host_result_Ellipses->dev_available_spot,max_size_global_result*sizeof(bool),cudaMemcpyDeviceToHost);
	else return host_width*host_height;

	unsigned int count = 0;
	for (unsigned int i=0;i<max_size_global_result;i++)
		if (!host_spots_available[i])   count ++;

	nb_Objects_found = count;
	return count;
}

extern "C" void getResult_Ellipse_Shape(Ellipse_Shape** a, int * b, double *c, bool verbose)
{
	cudaThreadSynchronize();//make sure last call of a kernel ended

	nb_Objects_found = getNbResult();
	if (nb_Objects_found == 0)
	{
		*b = 0; return;
	}
	cudaMemcpy(host_spots_available, host_result_Ellipses->dev_available_spot,max_size_global_result*sizeof(bool),cudaMemcpyDeviceToHost);


	float* x = (float*) malloc(max_size_global_result*sizeof(float));
	float* y = (float*) malloc(max_size_global_result*sizeof(float));
	float* angle = (float*) malloc(max_size_global_result*sizeof(float));
	float* axis0 = (float*) malloc(max_size_global_result*sizeof(float));
	float* axis1 = (float*) malloc(max_size_global_result*sizeof(float));

	cudaMemcpy(x, host_result_Ellipses->dev_x,					max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);
	cudaMemcpy(y, host_result_Ellipses->dev_y,					max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);
	cudaMemcpy(angle, host_result_Ellipses->dev_angle,			max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);
	cudaMemcpy(axis0, host_result_Ellipses->dev_axis0,			max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);
	cudaMemcpy(axis1, host_result_Ellipses->dev_axis1,			max_size_global_result*sizeof(float),cudaMemcpyDeviceToHost);

	float localEnergy;
	unsigned int* Pixel_i = (unsigned int*) malloc((host_width+2*host_border)*(host_height+2*host_border)*sizeof(unsigned int));
	unsigned int* minArea			  = (unsigned int*) malloc((host_width+2*host_border)*(host_height+2*host_border)*sizeof(unsigned int)*NB_OVERLAP_MAX_ON_PIXEL);
	
	cudaMemcpy(Pixel_i, dev_nbObjectsAtPixel_i,(host_width+2*host_border)*(host_height+2*host_border)*sizeof(unsigned int)     ,cudaMemcpyDeviceToHost);
	cudaMemcpy(minArea,	            dev_minArea,   		    (host_width+2*host_border)*(host_height+2*host_border)*sizeof(unsigned int)*NB_OVERLAP_MAX_ON_PIXEL,cudaMemcpyDeviceToHost);


	*a = (Ellipse_Shape*) malloc(nb_Objects_found*sizeof(Ellipse_Shape));
	unsigned int count = 0;
	double finalEnergy = 0;
	for (unsigned int i=0;i<max_size_global_result/* && count < nb_Objects_found*/;i++)
	{
		if (!host_spots_available[i])
		{

			localEnergy = getEnergy(x[i],y[i],angle[i],axis0[i],axis1[i],Pixel_i,minArea, verbose, 0);
			(*a)[count++] = Ellipse_Shape(x[i],y[i],angle[i],axis0[i],axis1[i]);
			finalEnergy += localEnergy;

			if (host_param->Verbose > 1)printf("%d energy %f %f\n",count,localEnergy, finalEnergy);
		}
	}

	*b = count;//nb_Objects_found;
	*c = finalEnergy;

	free(x);free(y);free(angle);free(axis0);free(axis1);

	free(Pixel_i);free(minArea);
}

extern "C" void cuda_startTimer()
{
	cudaEventCreate(&start_event);
	cudaEventCreate(&stop_event);
	cudaEventRecord(start_event, 0);
}

extern "C" void cuda_resetTime()
{
	cudaEventRecord(start_event, 0);
}

extern "C" void cuda_blockandWait()
{
	cudaThreadSynchronize();
}

extern "C" float cuda_getTime()
{
	//record the time 
	cudaEventRecord(stop_event, 0);
	cudaEventSynchronize(stop_event);
	float time;
	cudaEventElapsedTime( &time, start_event, stop_event );
	cuda_resetTime();
	return time;
}

extern "C" void cuda_delete()
{

	if (dev_rand) cudaFree(dev_rand);
	if (dev_BirthMap) cudaFree(dev_BirthMap);
	if (dev_dataAtPixel)  cudaFree(dev_dataAtPixel);
	if (dev_minArea) cudaFree(dev_minArea);
	if (dev_nbObjectsAtPixel_i) cudaFree(dev_nbObjectsAtPixel_i);
	if (dev_nbObjectsAtCell_i) cudaFree(dev_nbObjectsAtCell_i);

	if (dev_image) 
		cudaFree(dev_image);

	if (dev_param) cudaFree(dev_param);
   	if (host_result_Ellipses)
	{
		if (host_result_Ellipses->dev_x)cudaFree(host_result_Ellipses->dev_x);
		if (host_result_Ellipses->dev_y)cudaFree(host_result_Ellipses->dev_y);
		if (host_result_Ellipses->dev_angle)cudaFree(host_result_Ellipses->dev_angle);
		if (host_result_Ellipses->dev_axis0)cudaFree(host_result_Ellipses->dev_axis0);
		if (host_result_Ellipses->dev_axis1)cudaFree(host_result_Ellipses->dev_axis1);
		if (host_result_Ellipses->dev_available_spot)cudaFree(host_result_Ellipses->dev_available_spot);
		if (host_result_Ellipses->dev_index_available_spots)cudaFree(host_result_Ellipses->dev_index_available_spots);
		if (host_result_Ellipses->dev_index_not_available_spots)cudaFree(host_result_Ellipses->dev_index_not_available_spots);
		free(host_result_Ellipses);
	}
	if (dev_result_Ellipses)cudaFree(dev_result_Ellipses);
	if (host_index_spots_available)	free(host_index_spots_available);
	if (host_spots_available) free(host_spots_available);

	isInit = false;

	if (host_quadtree)
	{
		delete host_quadtree;   host_quadtree = NULL;
		if (dev_array_quadtree_x) cudaFree(dev_array_quadtree_x);
		if (dev_array_quadtree_y) cudaFree(dev_array_quadtree_y);
		if (dev_array_quadtree_w) cudaFree(dev_array_quadtree_w);
		if (dev_array_quadtree_h) cudaFree(dev_array_quadtree_h);
	}


	if (manifoldTexture){ free(manifoldTexture);manifoldTexture = NULL;}; 
	host_quadtree = NULL;

	curandDestroyGenerator(gen);
}

void init_Result_Ellipse(int size);

unsigned int initGPU_RJMCMC_image(int* ptrImage, int w, int h, int tmp_size_cell,Const_Parameter* param, int border)
{
	host_image = ptrImage;
	host_param = param;
	host_width = w;
	host_height = h;
	host_border = border;
	size_cell = tmp_size_cell; 


	//---------------------just for display purpose-----------------------------
	int iddevice = 0;//1 quadro, 0 GTX 480			//<----------- select GPU here
	unsigned int free = initDevice(iddevice);
	//---------------------------------------------------------------------------


	cudaMalloc(&dev_minArea,(w+2*border)*(h+2*border)*sizeof(unsigned int)*NB_OVERLAP_MAX_ON_PIXEL);	

	cudaMalloc(&dev_nbObjectsAtPixel_i,	(w+2*border)*(h+2*border)*sizeof(unsigned int));
	cudaMemset(dev_nbObjectsAtPixel_i,0,	(w+2*border)*(h+2*border)*sizeof(unsigned int));
		
	int nb_rnd_val = 6;
	dev_param = NULL;
	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{
	cudaMalloc(&dev_param,					sizeof(bayesian_Parameter));
	cudaMemcpy(dev_param, param,			sizeof(bayesian_Parameter), cudaMemcpyHostToDevice);
	}

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA )
	{
		printf("in init	for batha (ellipse)\n");
	cudaMalloc(&dev_param,					sizeof(bhatt_Parameter));
	cudaMemcpy(dev_param, param,			sizeof(bhatt_Parameter), cudaMemcpyHostToDevice);
	}

	for (int i=0;i<4;i++)
	{
		nb_cell_grid[i][0] = ceil( 0.5*(float(w+2*border)/size_cell ) ) ; 
		nb_cell_grid[i][1] = ceil( 0.5*(float(h+2*border)/size_cell ) ) ;
		nb_Cells = nb_cell_grid[i][0]*nb_cell_grid[i][1];
	}

	cudaMalloc(&dev_nbObjectsAtCell_i,		4*nb_Cells*sizeof(unsigned int));
	cudaMemset(dev_nbObjectsAtCell_i,0,	4*nb_Cells*sizeof(unsigned int));

	nb_Objects_found = 0;
	max_size_global_result = X_INCREASE_MEMORY*4*nb_Cells;														    



	init_Result_Ellipse(max_size_global_result);

	//Image on dev
	dev_image = NULL;
	cudaMalloc(&dev_image,					w*h*sizeof(int));
	cudaMemcpy(dev_image, ptrImage,			w*h*sizeof(int), cudaMemcpyHostToDevice);


	//rand generator
	printf("init rand generator : %s\n",curandCreateGenerator(&gen,CURAND_RNG_PSEUDO_DEFAULT) == CURAND_STATUS_SUCCESS ? "Done" : "failed");
	if (curandSetPseudoRandomGeneratorSeed(gen, 1234ULL) != CURAND_STATUS_SUCCESS){checkCUDAError("error in init");printf("(((EOR))) init rand generator\n");exit(1);};
	cudaMalloc(&dev_rand, nb_rnd_val*nb_Cells*sizeof(float));
	//warmup ? -- I read somewhere that the first generation take more time, so I call it here for 'warmup'
 	curandStatus_t status = curandGenerateUniform(gen, dev_rand, nb_Cells*nb_rnd_val);
	if (status != CURAND_STATUS_SUCCESS){checkCUDAError("error in init");printf("(((ERROR))) init rand generator %d\n",status);exit(1);}

	host_index_spots_available = (unsigned int*) malloc(nb_Cells*sizeof(unsigned int));
	host_spots_available = (bool*) malloc(max_size_global_result*sizeof(bool));

	return free;
}

//can be for ellipse or segment
extern "C" float cuda_RJMCMC_init(int* ptrImage, int w, int h, int half_size_cell, Const_Parameter* param,int border, METHOD_DATATERM method)  //0 ellipse, 1 segment
{
	isInit = true;
	method_dterm = method;

	unsigned int free_m = initGPU_RJMCMC_image(ptrImage,w,h,2*half_size_cell,param,border);
	printf("image loaded %d x %d, nbCells %d size_cell %d\n",w, h, nb_Cells,2*half_size_cell);

	ratio_birth_death = 0.5;
	if (param->ConstructDataMap)
	{	
		if (strlen(param->FileName) > 0)
		{
			bool force_square = true;
			cimg_library::CImg<int> inputImage(param->FileName); 
			printf("create quadtree, noHo !\n");

			host_quadtree = new quadtree(static_cast<int*>(inputImage.data()), inputImage.width(), inputImage.height(), 2*half_size_cell, false, border, force_square);
			host_quadtree->prepare();

			//Only for display purpose
			//host_quadtree->draw_pixelMap_pretty("data/", "sample_big_bad_quadtree");

			ratio_birth_death = host_quadtree->get_avg_cell();
			ratio_birth_death = 1-(ratio_birth_death/ (ratio_birth_death+nb_Cells));

			if (force_square)
			{
			   	cudaMalloc( &dev_array_quadtree_x, host_quadtree->arrayx.size()*sizeof(unsigned int) );
				cudaMemcpy(dev_array_quadtree_x, &(host_quadtree->arrayx[0]),host_quadtree->arrayx.size()*sizeof(unsigned int),cudaMemcpyHostToDevice);
				dev_array_quadtree_y = NULL;
				dev_array_quadtree_w = NULL;
				dev_array_quadtree_h = NULL;

			}else{
				cudaMalloc( &dev_array_quadtree_x, host_quadtree->arrayx.size()*sizeof(unsigned int) );
				cudaMemcpy(dev_array_quadtree_x, &(host_quadtree->arrayx[0]),host_quadtree->arrayx.size()*sizeof(unsigned int),cudaMemcpyHostToDevice);

				cudaMalloc( &dev_array_quadtree_y, host_quadtree->arrayy.size()*sizeof(unsigned int) );
				cudaMemcpy(dev_array_quadtree_y, &(host_quadtree->arrayy[0]),host_quadtree->arrayy.size()*sizeof(unsigned int),cudaMemcpyHostToDevice);

				cudaMalloc( &dev_array_quadtree_w, host_quadtree->arrayw.size()*sizeof(unsigned int) );
				cudaMemcpy(dev_array_quadtree_w, &(host_quadtree->arrayw[0]),host_quadtree->arrayw.size()*sizeof(unsigned int),cudaMemcpyHostToDevice);

				cudaMalloc( &dev_array_quadtree_h, host_quadtree->arrayh.size()*sizeof(unsigned int) );
				cudaMemcpy(dev_array_quadtree_h, &(host_quadtree->arrayh[0]),host_quadtree->arrayh.size()*sizeof(unsigned int),cudaMemcpyHostToDevice);
			}

		}
	}

	unsigned int free2,total;
	cuMemGetInfo((size_t*)&free2,(size_t*) &total);
	printf("memory left: %f / %f (Mb)\n",(total - free2)/ 1024. / 1024., total/ 1024. / 1024.);
	checkCUDAError("error at init");
	printf("end init cuda RJMCMC (#Threads lower level %d) (ratio b/d %f)\n", nb_Cells,ratio_birth_death);
	return ratio_birth_death;
}


//----------------------------------------------------------------------------

//------------------------------------------------Ellipse
void init_Result_Ellipse(int size)
{
	host_result_Ellipses = (Ellipse_result*) malloc (sizeof(Ellipse_result));

	cudaMalloc(&(host_result_Ellipses->dev_x)							, size*sizeof(float));
	cudaMalloc(&(host_result_Ellipses->dev_y)							, size*sizeof(float));
	cudaMalloc(&(host_result_Ellipses->dev_angle)						, size*sizeof(float));
	cudaMalloc(&(host_result_Ellipses->dev_axis0)						, size*sizeof(float));
	cudaMalloc(&(host_result_Ellipses->dev_axis1)						, size*sizeof(float));
	cudaMalloc(&(host_result_Ellipses->dev_available_spot)				, size*sizeof(bool));

	{
		cudaMalloc(&(host_result_Ellipses->dev_index_available_spots)		, nb_Cells*sizeof(unsigned int));

		unsigned int* tmp = (unsigned int*) malloc(nb_Cells*sizeof(unsigned int));
		for (int i=0;i<nb_Cells;i++)
			tmp[i] = i;
		cudaMemcpy( host_result_Ellipses->dev_index_available_spots , tmp , nb_Cells*sizeof(unsigned int), cudaMemcpyHostToDevice );
		free(tmp);

	}

	cudaMalloc(&(host_result_Ellipses->dev_index_not_available_spots)	, size*sizeof(int));//use for death

	cudaMemset(host_result_Ellipses->dev_available_spot,true			, size*sizeof(bool));

	cudaMalloc( (void **)&(dev_result_Ellipses) , sizeof(Ellipse_result));
	cudaMemcpy( dev_result_Ellipses , host_result_Ellipses, sizeof(Ellipse_result), cudaMemcpyHostToDevice );


}
void increaseMemResult_Ellipse()
{
	//NOTE (TODO): That function is bugged, if the code go in, it will crash later (event if it exits without error this function) , so need to avoid calling it for now
	cudaThreadSynchronize();
	if (host_param->Verbose > 1) printf("in increase\n");

	cudaFree(dev_result_Ellipses);
	Ellipse_result *old_p_host = host_result_Ellipses;

	init_Result_Ellipse(max_size_global_result+X_INCREASE_MEMORY*nb_Cells);//reinit dev and host bigger

	cudaMemcpy(host_result_Ellipses->dev_x, old_p_host->dev_x,													max_size_global_result*sizeof(float)			,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_y, old_p_host->dev_y,													max_size_global_result*sizeof(float)			,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_angle, old_p_host->dev_angle,											max_size_global_result*sizeof(float)			,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_axis0, old_p_host->dev_axis0,											max_size_global_result*sizeof(float)			,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_axis1, old_p_host->dev_axis1,											max_size_global_result*sizeof(float)			,cudaMemcpyDeviceToDevice);

	cudaMemcpy(host_result_Ellipses->dev_available_spot, old_p_host->dev_available_spot,						max_size_global_result*sizeof(bool)				,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_index_not_available_spots, old_p_host->dev_index_not_available_spots,	max_size_global_result*sizeof(int)				,cudaMemcpyDeviceToDevice);
	cudaMemcpy(host_result_Ellipses->dev_index_available_spots, old_p_host->dev_index_available_spots,			nb_Cells*sizeof(unsigned int)					,cudaMemcpyDeviceToDevice);


	cudaFree(old_p_host->dev_x);cudaFree(old_p_host->dev_y);cudaFree(old_p_host->dev_angle);
	cudaFree(old_p_host->dev_axis0);cudaFree(old_p_host->dev_axis1);
	cudaFree(old_p_host->dev_available_spot);cudaFree(old_p_host->dev_index_available_spots);
	cudaFree(old_p_host->dev_index_not_available_spots);
	free(old_p_host);	

	free(host_spots_available);
	host_spots_available = (bool*) malloc(max_size_global_result*sizeof(bool));

	max_size_global_result+=X_INCREASE_MEMORY*nb_Cells;

	if (host_param->Verbose > 1) printf("out increase\n");

}
void updateAvailableSpot_Ellipse()	 
{
	if (host_param->Verbose > 1) printf("in update Ellipse\n");
	cudaMemcpy(host_spots_available, host_result_Ellipses->dev_available_spot,max_size_global_result*sizeof(bool),cudaMemcpyDeviceToHost);

	int count_available = 0;
	int count_not_available = 0;
	unsigned int last_index = max_size_global_result;


	for (unsigned int i=0;i<max_size_global_result;i++)
	{
		if (host_spots_available[i])
		{
			if (count_available < nb_Cells)
				host_index_spots_available[count_available++] = i;
		}else{
			count_not_available++;
			last_index = i;
		}
	}


	//complete the array for available
	for (int i=0;i<nb_Cells-count_available;i++)
		host_index_spots_available[i+count_available] = i+last_index;

	nb_Objects_found = count_not_available;

	cudaMemcpy(host_result_Ellipses->dev_index_available_spots,     host_index_spots_available,     nb_Cells*sizeof(unsigned int), cudaMemcpyHostToDevice);


	if (host_index_spots_available[nb_Cells - 1] > max_size_global_result)//if the last index marked as available is bigger than the memory available, increase memory
		increaseMemResult_Ellipse();

}
extern "C" void birth_Simulation_GPU_Ellipse(float temperature, int* ptrImage, int w, int h, int border, Const_Parameter* param)
{
	if (!isInit){std::cerr << "must be init before "<<std::endl;return;}

	int grid2run; 
	int max_nb_threads;
	int nb_block;
	int level;

	if (host_param->Verbose > 1) printf("in birth (ellipse)\n");

	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{

		grid2run =  rand()%4;//in [0,3]

		max_nb_threads = nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];

		if ((status = curandGenerateUniform(gen, dev_rand,max_nb_threads*6)) != CURAND_STATUS_SUCCESS){printf("ABirth-> Could not set the number of outputs from number generator: %d\n",status);exit(1);}

		nb_block =  max(( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block,1);

		dim3 threads_size(host_nb_thread_per_block,1);
		dim3 blocks_size(nb_block,1 );

		Birth_on_Grid_RJMCMC_bayesian<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border, nb_Cells, //input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i, dev_nbObjectsAtCell_i,																	  //output
			temperature, (bayesian_Parameter*)dev_param, nb_Objects_found);

	}


	if (method_dterm == METHOD_DATATERM::BHATTACHARYA)
	{
		if (host_quadtree)
		{
			const unsigned int offset = host_quadtree->getOffset(0);	//0 default, 1 equiprobable
			const unsigned int nb_thread_total = host_quadtree->getNbCells(offset);
			const unsigned int positionStart = host_quadtree->getStartIndex(offset);
			const unsigned int size_cell_qt =  host_quadtree->getSizeCell(offset);
			level = offset/4;
			grid2run = offset - level*4;

			if ((status = curandGenerateUniform(gen, dev_rand,nb_thread_total*6)) != CURAND_STATUS_SUCCESS){printf("BBBirth-> Could not set the number of outputs from number generator: %d\n",status);exit(1);}

			nb_block =  max(( ( nb_thread_total  )+host_nb_thread_per_block-1)/host_nb_thread_per_block,1);

			dim3 threads_size(host_nb_thread_per_block,1);
			dim3 blocks_size(nb_block,1 );

			if (host_quadtree->isConstrained())
			{
				Birth_on_Grid_RJMCMC_bhat_nonHo<<<blocks_size, threads_size>>>
					(dev_image, w+2*border, h+2*border ,dev_rand, size_cell, grid[grid2run][0], grid[grid2run][1] , border, nb_Cells, //input
					dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i, dev_nbObjectsAtCell_i,	
					dev_array_quadtree_x, size_cell_qt, nb_thread_total, positionStart, ratio_birth_death,   //we put the grid index in dev_array_quadtree_x
					temperature, (bhatt_Parameter*) dev_param, nb_Objects_found);
			}else{
				printf("e am here, should I ?\n");
				Birth_on_Grid_RJMCMC_bhat_nonHo<<<blocks_size, threads_size>>>
					(dev_image, w+2*border, h+2*border ,dev_rand, size_cell, grid[grid2run][0], grid[grid2run][1] , border, nb_Cells, //input
					dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i, dev_nbObjectsAtCell_i,	
					dev_array_quadtree_x, dev_array_quadtree_y, dev_array_quadtree_w, dev_array_quadtree_h,  nb_thread_total, positionStart,
					temperature, (bhatt_Parameter*) dev_param, nb_Objects_found);
			}


			

		}
		else //no quadtree
		{
			grid2run = rand()%4;//in [0,3]

			max_nb_threads = nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];//nb_Cells;//nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];

			if ((status = curandGenerateUniform(gen, dev_rand,nb_Cells*6)) != CURAND_STATUS_SUCCESS){printf("error:: %s\n",cudaGetErrorString(cudaGetLastError()));printf("Birth-> Could not set the number of outputs from number generator: %d (%d %d %d)\n",status, &gen, dev_rand,max_nb_threads);exit(1);}
			nb_block =  max(( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block,1);

			dim3 threads_size(host_nb_thread_per_block,1);
			dim3 blocks_size(nb_block,1 );



			Birth_on_Grid_RJMCMC_bhat<<<blocks_size, threads_size>>>
				(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border, nb_Cells, //input
				dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i, dev_nbObjectsAtCell_i,																	  //output
				temperature, (bhatt_Parameter*) dev_param, nb_Objects_found);
		}
	}

	updateAvailableSpot_Ellipse();
};
extern "C" void death_Simulation_GPU_Ellipse(float temperature, int* ptrImage, int w, int h, int border, Const_Parameter* param)
{
	if (!isInit){std::cerr << "must be init before "<<std::endl;return;}

	if (nb_Objects_found == 0) return;
	if (host_param->Verbose > 1) printf("in death ellipse\n");
	int grid2run;
	int max_nb_threads;
	int nb_block;

	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{

		grid2run = rand()%4;//in [0,3]
		prepareEllipsesCandidateForGrid_i(grid2run,w,h);

		max_nb_threads = nb_Cells;
		if (curandGenerateUniform(gen, dev_rand, max_nb_threads) != CURAND_STATUS_SUCCESS){printf("Death: -> Could not set the number of outputs from number generator: %d\n", nb_Objects_found);exit(1);}


		nb_block =  max(( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block,1);

		dim3 threads_size(host_nb_thread_per_block,1);
		dim3 blocks_size(nb_block,1 );


		Death_on_Grid_RJMCMC<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
			temperature,(bayesian_Parameter*)  dev_param, nb_Objects_found, nb_Cells);

	}

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA)
	{
		grid2run = rand()%4;//in [0,3]
		prepareEllipsesCandidateForGrid_i(grid2run,w,h);

		max_nb_threads = nb_Cells;

		if (curandGenerateUniform(gen, dev_rand, max_nb_threads) != CURAND_STATUS_SUCCESS){printf("Death: -> Could not set the number of outputs from number generator: %d (%d %d %d)\n", nb_Objects_found, &gen, dev_rand, max_nb_threads);exit(1);}

		nb_block =  max(( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block,1);

		dim3 threads_size(host_nb_thread_per_block,1);
		dim3 blocks_size(nb_block,1 );


		Death_on_Grid_RJMCMC_bhat<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 0.5,																				//output
			 temperature,(bhatt_Parameter*)  dev_param, nb_Objects_found, nb_Cells);
		
	}

	updateAvailableSpot_Ellipse();
};
extern "C" void Translate_Simulation_GPU_Ellipse(float temperature, int* ptrImage, int w, int h, int border, Const_Parameter* param)
{
	if (!isInit){std::cerr << "must be init before "<<std::endl;return;}

	if (nb_Objects_found == 0) return;
	if (host_param->Verbose > 1) printf("in translate ellipse\n");
	int grid2run = rand()%4;//in [0,3]

	prepareEllipsesCandidateForGrid_i(grid2run,w,h);

	int max_nb_threads = nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];

	if ((status = curandGenerateUniform(gen, dev_rand,max_nb_threads*3)) != CURAND_STATUS_SUCCESS){printf("Birth-> Could not set the number of outputs from number generator: %d\n",status);exit(1);}

	int nb_block =  ( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block;

	dim3 threads_size(host_nb_thread_per_block,1);
	dim3 blocks_size(nb_block,1 );

	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{
		Translate_on_Grid<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
			 temperature,(bayesian_Parameter*) dev_param, nb_Objects_found, nb_Cells);																			//const 
	}

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA)
	{
			Translate_on_Grid_bhat<<<blocks_size, threads_size>>>
				(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
				dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
				temperature,(bhatt_Parameter*) dev_param, nb_Objects_found, nb_Cells);	
	}

	updateAvailableSpot_Ellipse();

};
extern "C" void Enlarge_Simulation_GPU_Ellipse(float temperature, int* ptrImage, int w, int h, int border, Const_Parameter* param)
{
	if (!isInit){std::cerr << "must be init before "<<std::endl;return;}

	if (nb_Objects_found == 0) return;
	if (host_param->Verbose > 1) printf("in enlarge ellipse\n");
	int grid2run = rand()%4;//in [0,3]

	prepareEllipsesCandidateForGrid_i(grid2run,w,h);

	int max_nb_threads = nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];


	if ((status = curandGenerateUniform(gen, dev_rand,max_nb_threads*3)) != CURAND_STATUS_SUCCESS){printf("Birth-> Could not set the number of outputs from number generator: %d\n",status);exit(1);}

	int nb_block =  ( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block;

	dim3 threads_size(host_nb_thread_per_block,1);
	dim3 blocks_size(nb_block,1 );

	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{
		EnlargeShrink_on_Grid<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
			 temperature,(bayesian_Parameter*) dev_param, nb_Objects_found, nb_Cells);																			//const 
	}

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA)
	{

			EnlargeShrink_on_Grid_bhat<<<blocks_size, threads_size>>>
				(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
				dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
				 temperature,(bhatt_Parameter*) dev_param, nb_Objects_found, nb_Cells);
	}
};
extern "C" void Rotate_Simulation_GPU_Ellipse(float temperature, int* ptrImage, int w, int h, int border, Const_Parameter* param)
{
	if (!isInit){std::cerr << "must be init before "<<std::endl;return;}

	if (nb_Objects_found == 0) return;
	if (host_param->Verbose > 1) printf("in rotate ellipse\n");
	int grid2run = rand()%4;//in [0,3]

	prepareEllipsesCandidateForGrid_i(grid2run,w,h);

	int max_nb_threads = nb_cell_grid[grid2run][0]*nb_cell_grid[grid2run][1];


	if ((status = curandGenerateUniform(gen, dev_rand,max_nb_threads*2)) != CURAND_STATUS_SUCCESS){printf("Birth-> Could not set the number of outputs from number generator: %d\n",status);exit(1);}

	int nb_block =  ( ( max_nb_threads  )+host_nb_thread_per_block-1)/host_nb_thread_per_block;

	dim3 threads_size(host_nb_thread_per_block,1);
	dim3 blocks_size(nb_block,1 );

	if (method_dterm == METHOD_DATATERM::BAYESIAN)
	{
		Rotate_on_Grid<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
			 temperature,(bayesian_Parameter*) dev_param, nb_Objects_found, nb_Cells);
	} 

	if (method_dterm == METHOD_DATATERM::BHATTACHARYA)
	{
		    Rotate_on_Grid_bhat<<<blocks_size, threads_size>>>
			(dev_image, w+2*border, h+2*border ,dev_rand, nb_cell_grid[grid2run][0], nb_cell_grid[grid2run][1], size_cell,  grid[grid2run][0], grid[grid2run][1], border,  //input																				//input
			dev_result_Ellipses, dev_minArea, dev_nbObjectsAtPixel_i,dev_nbObjectsAtCell_i, 																				//output
			 temperature,(bhatt_Parameter*) dev_param, nb_Objects_found, nb_Cells);
	} 

};





///////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////
















