 
 		*ret_index = best_idx;
 		*ret_dist = best_dist;
+	}
+
+	__device__ void SearchAtNode_knn(const CUDA_KDNode *nodes, const int *indexes, const Kd_tree_point *pts, int cur, const Kd_tree_point &query, int *ret_index, float *ret_dist, int *ret_node, int k)
+	{
+		int neighbor_nb = 0;
+
+
+		//travasel to the nodes
+		while (true)
+		{
+			int split_axis = nodes[cur].level % KDTREE_DIM;
+
+			if (nodes[cur].left == -1){
+				//Get to the leaf node
+
+				neighbor_nb += nodes[cur].num_indexes;
+				while (neighbor_nb < k)
+				{
+					cur = nodes[cur].parent;
+					neighbor_nb = nodes[cur].num_indexes;
+				}
+				*ret_node = cur;
+
+				//Now we get enough neighbors in cur node
+				int *temp_index = (int*)malloc(sizeof(int) * nodes[cur].num_indexes);
+				float *temp_dists = (float*)malloc(sizeof(float) * nodes[cur].num_indexes);
+
+				for (int i = 0; i < nodes[cur].num_indexes; ++i)
+				{
+					temp_index[i] = indexes[nodes[cur].indexes + i];
+					temp_dists[i] = Distance(query, pts[temp_index[i]]);
+				}
+
+				int n = nodes[cur].num_indexes;
+				//kðݵõǰСľ
+				int best_idx = 0;
+				double best_dist = FLT_MAX;
+				for (int i = 0; i < k; ++i)
+				{
+					for (int j = i; j < n; ++j)
+					{
+						if (temp_dists[j] < best_dist)
+						{
+							best_dist = temp_dists[j];
+							best_idx = temp_index[j];
+
+							temp_dists[j] = temp_dists[i];
+							temp_index[j] = temp_index[i];
+
+							temp_index[i] = best_idx;
+							temp_dists[i] = best_dist;
+						}
+					}
+
+					ret_index[i] = best_idx;
+					ret_dist[i] = best_dist;
+
+					best_idx = 0;
+					best_dist = FLT_MAX;
+				}
+				break;
+			}
+			else if (query.coords[split_axis] < nodes[cur].split_value){
+				cur = nodes[cur].left;
+			}
+			else{
+				cur = nodes[cur].right;
+			}
+		}
+
+	}
+
+	__device__ void SearchAtiNodeRange_knn(const CUDA_KDNode *nodes, const int *indexes, const Kd_tree_point *pts, const Kd_tree_point &query, int cur, double range, int* &ret_index, double* &ret_dist, int k)
+	{
+		// Goes through all the nodes that within "radius"
+		int to_visit[10];
+		int to_visit_pos = 0;
+
+		to_visit[to_visit_pos++] = cur;
+
+		while (to_visit_pos){
+			int next_search[10];
+			int next_search_pos = 0;
+
+			while (to_visit_pos){
+				cur = to_visit[to_visit_pos - 1];
+				to_visit_pos--;
+
+				int split_axis = nodes[cur].level % KDTREE_DIM;
+
+				//already get to the leaf
+				if (nodes[cur].left == -1){
+					for (int i = 0; i < nodes[cur].num_indexes; ++i)
+					{
+						int idx = indexes[nodes[cur].indexes + i];
+						float d = Distance(query, pts[idx]);
+
+						for (int j = 0; j < k; ++j)
+						{
+
+							if (d < ret_dist[j])
+							{
+								for (int q = k - 1; q > j; --q)
+								{
+									ret_index[q] = ret_index[q - 1];
+									ret_dist[q] = ret_dist[q - 1];
+								}
+								ret_index[j] = idx;
+								ret_dist[j] = d;
+								break;
+							}
+						}
+					}
+				}
+				else{
+					double d = query.coords[split_axis] - nodes[cur].split_value;
+
+					// There are 3 possible scenarios
+					// The hypercircle only intersects the left region
+					// The hypercircle only intersects the right region
+					// The hypercricle intersects both
+
+					if (fabs(d) > range) {
+						if (d < 0)
+							next_search[next_search_pos++] = nodes[cur].left;
+						else
+							next_search[next_search_pos++] = nodes[cur].right;
+					}
+					else {
+						next_search[next_search_pos++] = nodes[cur].left;
+						next_search[next_search_pos++] = nodes[cur].right;
+					}
+				}
+			}
+			// No memcpy available??
+			for (int i = 0; i < next_search_pos; i++)
+				to_visit[i] = next_search[i];
+
+			to_visit_pos = next_search_pos;
+		}
+	}
+
+
+	__device__ void Search_knn(const CUDA_KDNode *nodes, const int *indexes, const Kd_tree_point *pts, const Kd_tree_point &query, int *ret_index, double *ret_dist, int k)
+	{
+		// Find the first closest node, this will be the upper bound for the next searches
+		int best_node = 0;
+		int* k_idx = (int*)malloc(sizeof(int) * k);
+		double* k_dist = (double*)malloc(sizeof(float) * k);
+		double radius = 0;
+
+		SearchAtNode_knn(nodes, indexes, pts, 0 /* root */, query, k_idx, k_dist, &best_node, k);
+
+		radius = sqrt(k_dist[k - 1]);
+
+		//Now find other posiible candidates
+		int cur = best_node;
+
+		while (nodes[cur].parent != -1)
+		{
+			int parent = nodes[cur].parent;
+			int split_axis = nodes[cur].level % KDTREE_DIM;
+
+			//Search the other node
+
+			if (fabs(nodes[parent].split_value - query.coords[split_axis]) <= radius)
+			{
+
+				if (nodes[parent].left != cur)
+					SearchAtiNodeRange_knn(nodes, indexes, pts, query, nodes[parent].left, radius, k_idx, k_dist, k);
+				else
+					SearchAtiNodeRange_knn(nodes, indexes, pts, query, nodes[parent].right, radius, k_idx, k_dist, k);
+
+			}
+			cur = parent;
+		}
+
+		for (int i = 0; i < k; ++i)
+		{
+			ret_index[i] = k_idx[i];
+			ret_dist[i] = k_dist[i];
+		}
+		free(k_idx);
+		free(k_dist);
+	}
+
+
+	__global__ void SearchBatch_knn(const CUDA_KDNode *nodes, const int *indexes, const Kd_tree_point *pts, int num_pts, Kd_tree_point *queries, int num_queries, int *ret_index, float *ret_dist, int k)
+	{
+		int idx = blockIdx.x*blockDim.x + threadIdx.x;
+
+		if (idx >= num_queries)
+			return;
+
+		//test(nodes, indexes, pts, queries[idx], &ret_index[idx * k], &ret_dist[idx * k], k);
+		Search_knn(nodes, indexes, pts, queries[idx], &ret_index[idx * k], &ret_dist[idx * k], k);
 	}
 
 	__global__ void SearchBatch(const CUDA_KDNode *nodes, const int *indexes, const Kd_tree_point *pts, int num_pts, Kd_tree_point *queries, int num_queries, int *ret_index, double *ret_dist)

		cudaFree(gpu_ret_dist);

	}
}
		CheckCUDAError("Search");

		cudaMemcpy(&indexes[0], gpu_ret_indexes, sizeof(int)*queries.size(), cudaMemcpyDeviceToHost);
		cudaMemcpy(&dists[0], gpu_ret_dist, sizeof(double)*queries.size(), cudaMemcpyDeviceToHost);

		cudaFree(gpu_queries);
		cudaFree(gpu_ret_indexes);
