 		//Mat1i init_labels;
 		//d_block_labels_.download(init_labels);
 
-		Merge << <grid_size_, block_size_ >> > (d_img_, d_img_labels_);
+		Merge << <grid_size_, block_size_ >> > (d_img_labels_, last_pixel_);
 
 		//Mat1i block_info_final;
 		//d_img_labels_.download(block_info_final);		

		grid_size_ = dim3((((d_img_.cols + 1) / 2) + BLOCK_COLS - 1) / BLOCK_COLS, (((d_img_.rows + 1) / 2) + BLOCK_ROWS - 1) / BLOCK_ROWS, 1);
		block_size_ = dim3(BLOCK_COLS, BLOCK_ROWS, 1);

		InitLabeling << <grid_size_, block_size_ >> > (d_img_, d_img_labels_);

		Compression << <grid_size_, block_size_ >> > (d_img_labels_);


	}


	__global__ void InitLabeling(const cuda::PtrStepSzb img, cuda::PtrStepSzi labels, unsigned char *last_pixel) {
		unsigned row = (blockIdx.y * BLOCK_ROWS + threadIdx.y) * 2;
		unsigned col = (blockIdx.x * BLOCK_COLS + threadIdx.x) * 2;
		unsigned img_index = row * img.step + col;
		unsigned labels_index = row * (labels.step / labels.elem_size) + col;
		
		if (row < labels.rows && col < labels.cols) {

			unsigned P = 0;			

			char buffer[4];
			*(reinterpret_cast<int*>(buffer)) = 0;

			if (col + 1 < img.cols) {
				// This does not depend on endianness
				*(reinterpret_cast<int16_t*>(buffer)) = *(reinterpret_cast<int16_t*>(img.data + img_index));

				if (row + 1 < img.rows) {
					*(reinterpret_cast<int16_t*>(buffer + 2)) = *(reinterpret_cast<int16_t*>(img.data + img_index + img.step));
				}
			}
			else {
				buffer[0] = img.data[img_index];

				if (row + 1 < img.rows) {
					buffer[2] = img.data[img_index + img.step];
				}
			}

			if (buffer[0]) {
				P |= 0x777;
			}
			if (buffer[1]) {
				P |= (0x777 << 1);
			}
			if (buffer[2]) {
				P |= (0x777 << 4);
			}

			if (col == 0) {
				P &= 0xEEEE;
			}
			if (col + 1 >= img.cols) {
				P &= 0x3333;
			}
			else if (col + 2 >= img.cols) {
				P &= 0x7777;
			}

			if (row == 0) {
				P &= 0xFFF0;
			}
			if (row + 1 >= img.rows) {
				P &= 0xFF;
			}
			else if (row + 2 >= img.rows) {
				P &= 0xFFF;
			}

			// P is now ready to be used to find neighbour blocks (or it should be)
			// P value avoids range errors
							
			int father_negative_offset = 0;
				
			// to_link is a bitmask that tells which squares are to be merged with the current one, in reduction phase.
			// First bit is North, Second bit is North-East and Third bit is West.
			unsigned char to_link = 0;							
			
			// North-West square
			if (HasBit(P, 0) && img.data[img_index - img.step - 1]) {
				father_negative_offset = 2 * (labels.step / labels.elem_size) + 2;
			}

			// North square
			if ((HasBit(P, 1) && img.data[img_index - img.step]) || (HasBit(P, 2) && img.data[img_index + 1 - img.step])) {
				if (!father_negative_offset) {
					father_negative_offset = 2 * (labels.step / labels.elem_size);
				}
				else {
					SetBit(to_link, 0);
				}
			}

			// North-East square
			if (HasBit(P, 3) && img.data[img_index + 2 - img.step]) {
				if (!father_negative_offset) {
					father_negative_offset = 2 * (labels.step / labels.elem_size) - 2;
				}
				else {
					SetBit(to_link, 1);
				}
			}

			// West square
			if ((HasBit(P, 4) && img.data[img_index - 1]) || (HasBit(P, 8) && img.data[img_index + img.step - 1])) {
				if (!father_negative_offset) {
					father_negative_offset = + 2;
				}
				else {
					SetBit(to_link, 2);
				}
			}

			// Just to see where to_link is stored
			SetBit(to_link, 4);
			
			labels.data[labels_index] = labels_index - father_negative_offset;
			if (col + 1 < labels.cols) {
				last_pixel = reinterpret_cast<unsigned char *>(labels.data + labels_index + 1);
			}
			else if (row + 1 < labels.rows) {
				last_pixel = reinterpret_cast<unsigned char *>(labels.data + labels_index + labels.step / labels.elem_size);
			}
			*last_pixel = to_link;
		}
	}

	__global__ void Merge(cuda::PtrStepSzi labels, unsigned char *last_pixel) {

		unsigned row = (blockIdx.y * BLOCK_ROWS + threadIdx.y) * 2;
		unsigned col = (blockIdx.x * BLOCK_COLS + threadIdx.x) * 2;
		unsigned labels_index = row * (labels.step / labels.elem_size) + col;

		if (row < labels.rows && col < labels.cols) {

			if (col + 1 < labels.cols) {
				last_pixel = reinterpret_cast<unsigned char *>(labels.data + labels_index + 1);
			}
			else if (row + 1 < labels.rows) {
				last_pixel = reinterpret_cast<unsigned char *>(labels.data + labels_index + labels.step / labels.elem_size);
			}
			unsigned char to_link = *last_pixel;

			if (HasBit(to_link, 0)) {
				Union(labels.data, labels_index, labels_index - 2 * (labels.step / labels.elem_size));
			}
			if (HasBit(to_link, 1)) {
				Union(labels.data, labels_index, labels_index - 2 * (labels.step / labels.elem_size) + 2);
			}
			if (HasBit(to_link, 2)) {
				Union(labels.data, labels_index, labels_index - 2);
			}
		}
	}

	__global__ void Compression(cuda::PtrStepSzi labels) {
		unsigned row = (blockIdx.y * BLOCK_ROWS + threadIdx.y) * 2;
		unsigned col = (blockIdx.x * BLOCK_COLS + threadIdx.x) * 2;
		unsigned labels_index = row * (labels.step / labels.elem_size) + col;

		if (row < labels.rows && col < labels.cols) {
			unsigned label = labels.data[labels_index];
			if (label != labels_index) {
				labels.data[labels_index] = Find(labels.data, label);
			}
		}
	}

	__global__ void FinalLabeling(const cuda::PtrStepSzb img, cuda::PtrStepSzi labels) {

