 		} while (assumed != old);
 
 		return __longlong_as_double(old);
+	}
+
+	/*
+	 * \brief fetch the double type in texture
+	 */
+	static __inline__ __device__
+		double fetch_double(texture<int2, 1> t, index_t i){
+		int2 v = tex1Dfetch(t, i);
+		return __hiloint2double(v.y, v.x);
 	}
 
 	/*


	__host__
		void CudaRestrictedVoronoiDiagram::print_return_data(std::string filename) const{
		/*for (int i = 0; i < points_nb_; ++i)
		{
			if (fabs(host_ret_[i * 4 + 3]) >= 1e-12){
				host_ret_[i * 4 + 0] /= host_ret_[i * 4 + 3];
				host_ret_[i * 4 + 1] /= host_ret_[i * 4 + 3];
				host_ret_[i * 4 + 2] /= host_ret_[i * 4 + 3];
			}
		}*/
		index_t line_num = 4;
		std::ofstream f;
		f.open(filename);

		}*/
	}

	/*
	* \brief Intersects a polygon with a points.
	*/
	__device__
		void c_intersection_clip_facet_SR(
		CudaPolygon& current_polygon,
		index_t i,
		index_t k
		){
		CudaPolygon polygon_buffer;

		//load /memory[points] 3 times.
		double3 pi = {
			c_points[i * 3 + 0],
			c_points[i * 3 + 1],
			c_points[i * 3 + 2]
		};

		for (index_t t = 0; t < k; ++t){

			//load /memory[points_nn] k times.
			index_t j = c_points_nn[i * k + t];

			if (i != j){
				//load /memroy[points] k * 3 times.
				double3 pj = {
					c_points[j * 3 + 0],
					c_points[j * 3 + 1],
					c_points[j * 3 + 2]
				};

				double dij = distance2(pi, pj);
				double R2 = 0.0;

				for (index_t tt = 0; tt < current_polygon.vertex_nb; ++tt){
					double3 pk = { current_polygon.vertex[tt].x, current_polygon.vertex[tt].y, current_polygon.vertex[tt].z };
					double dik = distance2(pi, pk);
					R2 = max(R2, dik);
				}
				if (dij > 4.1 * R2){
					return;
				}
				clip_by_plane(current_polygon, polygon_buffer, pi, pj, j);
				swap_polygon(current_polygon, polygon_buffer);
			}
		}
	}


	__global__
		void c_kernel(
		index_t			vertex_nb,
		index_t			points_nb,
		index_t*		facets, index_t			facets_nb,
		index_t			k_p,
		index_t*		facets_nn, index_t			dim,
		double*			retdata
		){
		index_t tid = blockIdx.x * blockDim.x + threadIdx.x;
		if (tid >= facets_nb) return;

		if (tid >= 0 && tid < points_nb){
			g_seeds_information[tid * 4 + 0] = 0.0;
			g_seeds_information[tid * 4 + 1] = 0.0;
			g_seeds_information[tid * 4 + 2] = 0.0;
			g_seeds_information[tid * 4 + 3] = 0.0;
		}

		//load \memory[facet] 3 times.
		int3 facet_index = {
			facets[tid * dim + 0],
			facets[tid * dim + 1],
			facets[tid * dim + 2]
		};

		//load \memory[vertex] 9 times.
		double3 v1 = {
			c_vertex[facet_index.x * dim + 0],
			c_vertex[facet_index.x * dim + 1],
			c_vertex[facet_index.x * dim + 2]
		};
		double3 v2 = {
			c_vertex[facet_index.y * dim + 0],
			c_vertex[facet_index.y * dim + 1],
			c_vertex[facet_index.y * dim + 2]
		};
		double3 v3 = {
			c_vertex[facet_index.z * dim + 0],
			c_vertex[facet_index.z * dim + 1],
			c_vertex[facet_index.z * dim + 2]
		};

		CudaPolygon current_polygon;
		current_polygon.vertex_nb = 3;

		current_polygon.vertex[0].x = v1.x; current_polygon.vertex[0].y = v1.y; current_polygon.vertex[0].z = v1.z; current_polygon.vertex[0].w = 1.0;
		current_polygon.vertex[1].x = v2.x; current_polygon.vertex[1].y = v2.y; current_polygon.vertex[1].z = v2.z; current_polygon.vertex[1].w = 1.0;
		current_polygon.vertex[2].x = v3.x; current_polygon.vertex[2].y = v3.y; current_polygon.vertex[2].z = v3.z; current_polygon.vertex[2].w = 1.0;

		CudaPolygon current_store = current_polygon;
		//doesn't have the stack?
		index_t to_visit[CUDA_Stack_size];
		index_t to_visit_pos = 0;

		index_t has_visited[CUDA_Stack_size];
		index_t has_visited_nb = 0;
		bool has_visited_flag = false;

		//load \memory[facets_nn] 1 time.
		to_visit[to_visit_pos++] = facets_nn[tid];
		has_visited[has_visited_nb++] = to_visit[0];

		while (to_visit_pos){
			index_t current_seed = to_visit[to_visit_pos - 1];
			to_visit_pos--;

			c_intersection_clip_facet_SR(
				current_polygon,
				current_seed,
				k_p
				);

			//now we get the clipped polygon stored in "polygon", do something.
			action(
			current_polygon,
			current_seed
			);

			//Propagate to adjacent seeds
			for (index_t v = 0; v < current_polygon.vertex_nb; ++v)
			{
				CudaVertex ve = current_polygon.vertex[v];
				int ns = ve.neigh_s;
				if (ns != -1)
				{
					for (index_t ii = 0; ii < has_visited_nb; ++ii)
					{
						//if the neighbor seed has clipped the polygon
						//the flag should be set "true"
						if (has_visited[ii] == ns)
							has_visited_flag = true;
					}
					//the neighbor seed is new!
					if (!has_visited_flag)
					{
						to_visit[to_visit_pos++] = ns;
						has_visited[has_visited_nb++] = ns;
					}
					has_visited_flag = false;
				}
			}
			current_polygon = current_store;
		}
		__syncthreads();


		/*for (index_t i = 0; i < points_nb; ++i){
		retdata[i] = g_seeds_polygon_nb[i];
		}*/
		for (index_t i = 0; i < points_nb * 4; ++i){
		retdata[i] = g_seeds_information[i];
		}
	}

	__device__
		void t_intersection_clip_facet_SR(
		CudaPolygon& current_polygon,
		index_t i,
		index_t k
		){
		CudaPolygon polygon_buffer;
		
		//load /memory[points] 3 times.
		double3 pi = {
			fetch_double(t_points, (int)(i * 3 + 0)),
			fetch_double(t_points, i * 3 + 1),
			fetch_double(t_points, i * 3 + 2)
		};
		
		for (index_t t = 0; t < k; ++t){

			//load /memory[points_nn] k times.
			index_t j = tex1Dfetch(t_points_nn, i * k + t);

			if (i != j){
				//load /memroy[points] k * 3 times.
				double3 pj = {
					fetch_double(t_points, j * 3 + 0),
					fetch_double(t_points, j * 3 + 1),
					fetch_double(t_points, j * 3 + 2)
				};

				double dij = distance2(pi, pj);
				double R2 = 0.0;

				for (index_t tt = 0; tt < current_polygon.vertex_nb; ++tt){
					double3 pk = { current_polygon.vertex[tt].x, current_polygon.vertex[tt].y, current_polygon.vertex[tt].z };
					double dik = distance2(pi, pk);
					R2 = max(R2, dik);
				}
				if (dij > 4.1 * R2){
					return;
				}
				clip_by_plane(current_polygon, polygon_buffer, pi, pj, j);
				swap_polygon(current_polygon, polygon_buffer);
			}
		}
	}

	__global__
		void t_kernel(
		index_t			vertex_nb,
		index_t			points_nb,
		index_t*		facets, index_t			facets_nb,
		index_t			k_p,
		index_t*		facets_nn, index_t			dim,
		double*			retdata
		){
		index_t tid = blockIdx.x * blockDim.x + threadIdx.x;
		if (tid >= facets_nb) return;

		if (tid >= 0 && tid < points_nb){
			g_seeds_information[tid * 4 + 0] = 0.0;
			g_seeds_information[tid * 4 + 1] = 0.0;
			g_seeds_information[tid * 4 + 2] = 0.0;
			g_seeds_information[tid * 4 + 3] = 0.0;
		}

		//load \memory[facet] 3 times.
		int3 facet_index = {
			facets[tid * dim + 0],
			facets[tid * dim + 1],
			facets[tid * dim + 2]
		};

		//load \memory[vertex] 9 times.
		double3 v1 = {
			fetch_double(t_vertex, facet_index.x * dim + 0),
			fetch_double(t_vertex, facet_index.x * dim + 1),
			fetch_double(t_vertex, facet_index.x * dim + 2)
		};
		double3 v2 = {
			fetch_double(t_vertex, facet_index.y * dim + 0),
			fetch_double(t_vertex, facet_index.y * dim + 1),
			fetch_double(t_vertex, facet_index.y * dim + 2)
		};
		double3 v3 = {
			fetch_double(t_vertex, facet_index.z * dim + 0),
			fetch_double(t_vertex, facet_index.z * dim + 1),
			fetch_double(t_vertex, facet_index.z * dim + 2)
		};

		CudaPolygon current_polygon;
		current_polygon.vertex_nb = 3;

		current_polygon.vertex[0].x = v1.x; current_polygon.vertex[0].y = v1.y; current_polygon.vertex[0].z = v1.z; current_polygon.vertex[0].w = 1.0;
		current_polygon.vertex[1].x = v2.x; current_polygon.vertex[1].y = v2.y; current_polygon.vertex[1].z = v2.z; current_polygon.vertex[1].w = 1.0;
		current_polygon.vertex[2].x = v3.x; current_polygon.vertex[2].y = v3.y; current_polygon.vertex[2].z = v3.z; current_polygon.vertex[2].w = 1.0;

		CudaPolygon current_store = current_polygon;
		//doesn't have the stack?
		index_t to_visit[CUDA_Stack_size];
		index_t to_visit_pos = 0;

		index_t has_visited[CUDA_Stack_size];
		index_t has_visited_nb = 0;
		bool has_visited_flag = false;

		//load \memory[facets_nn] 1 time.
		to_visit[to_visit_pos++] = facets_nn[tid];
		has_visited[has_visited_nb++] = to_visit[0];

		while (to_visit_pos){
			index_t current_seed = to_visit[to_visit_pos - 1];
			to_visit_pos--;

			t_intersection_clip_facet_SR(
				current_polygon,
				current_seed,
				k_p
				);

			//now we get the clipped polygon stored in "polygon", do something.
			/*action(
				current_polygon,
				current_seed
				);*/

			//Propagate to adjacent seeds
			for (index_t v = 0; v < current_polygon.vertex_nb; ++v)
			{
				CudaVertex ve = current_polygon.vertex[v];
				int ns = ve.neigh_s;
				if (ns != -1)
				{
					for (index_t ii = 0; ii < has_visited_nb; ++ii)
					{
						//if the neighbor seed has clipped the polygon
						//the flag should be set "true"
						if (has_visited[ii] == ns)
							has_visited_flag = true;
					}
					//the neighbor seed is new!
					if (!has_visited_flag)
					{
						to_visit[to_visit_pos++] = ns;
						has_visited[has_visited_nb++] = ns;
					}
					has_visited_flag = false;
				}
			}
			current_polygon = current_store;
		}
		__syncthreads();

		/*for (index_t i = 0; i < points_nb; ++i){
		retdata[i] = g_seeds_polygon_nb[i];
		}*/
		/*for (index_t i = 0; i < points_nb * 4; ++i){
			retdata[i] = g_seeds_information[i];
		}*/
	}

	__host__
	void CudaRestrictedVoronoiDiagram::compute_Rvd(){
		CudaStopWatcher watcher("compute_rvd");
		watcher.start();

		mode_ = TEXTURE_MEMORY;
		allocate_and_copy(mode_);
		//might be improved dim3 type.
		int threads = 512;
		int blocks = facet_nb_ / threads + ((facet_nb_ % threads) ? 1 : 0);
		switch (mode_)
		{
		case GLOBAL_MEMORY:
			kernel << < threads, blocks >> > (
				dev_vertex_, vertex_nb_,
				dev_points_, points_nb_,
				dev_facets_, facet_nb_,
				dev_points_nn_, k_,
				dev_facets_nn_, dimension_,
				dev_ret_
				);
			break;
		case CONSTANT_MEMORY:
			c_kernel << < threads, blocks >> >(
				vertex_nb_,
				points_nb_,
				dev_facets_, facet_nb_,
				k_,
				dev_facets_nn_, dimension_,
				dev_ret_
				);
			break;
		case TEXTURE_MEMORY:
			t_kernel << < threads, blocks >> >(
				vertex_nb_,
				points_nb_,
				dev_facets_, facet_nb_,
				k_,
				dev_facets_nn_, dimension_,
				dev_ret_
				);
			break;
		default:
			break;
		}

		CheckCUDAError("kernel function");
		copy_back();
		watcher.stop();
		watcher.synchronize();
		watcher.print_elaspsed_time(std::cout);
		
		//std::string out_file("..//out//S2_points_tex.txt");
		//print_return_data(out_file);
		free_memory();
	}

