__global__ void fill_lin_hough_array(float *ima_d,
						int imW, int imH,
						float *Sin_d, float *Cos_d,
						int Rmax, int angles, float treshold,
						float *hough_d){
	int xi = blockIdx.x * blockDim.x + threadIdx.x;
	int yi = blockIdx.y * blockDim.y + threadIdx.y;
	int i = xi + imW * yi;
	float x = (float)xi;
	float y = (float)yi;
	int Y = 0, k, R;
	if(xi >= imW || yi >= imH) return;
	float ima = ima_d[i];
	//float tresh = *treshold;
	if(ima > treshold){
		for(k = 0; k < angles; k++, Y+=Rmax){
			// R = x*cos(theta) + y*sin(theta)
			R = (int)(0.5f + x * Cos_d[k] + y * Sin_d[k]);
			if(R > 0 && R < Rmax) hough_d[R + Y] += ima;
		}
	}
}
/* 182мс
 * построение преобразования Хафа для поиска линий
 * вход:
 * 		ima - данные изображения
 * 		imW, imH - ширина и высота изображения
 * 		Rmax - максимальная величина R
 * 		angles - размер массива углов (шаг по углу равен 180/angles градусов)
 * выход: hough - инициализированный внешней функцией массив,
 * 		в котором будет размещено преобразование Хафа
 * 		!!! массив должен быть проинициализирован нулями ДО вызова этой функции
 * Выходной массив уже нормирован на единицу
 */
extern "C" int fill_hough_lines(float *ima, int imW, int imH, int Rmax, int angles, float *hough){
	#undef RETMACRO
	#define RETMACRO do{ ret = 0; goto free_all; }while(0)
	int sz, ret = 1;
	int lblksz = LBLKSZ;
	float *ima_d = NULL, *hough_d = NULL;
	float treshold = 0.5f;
	init_sincos(angles);
	sz = imW * imH;
	dim3 blkdim(QBLKSZ, QBLKSZ);
	dim3 griddim((imW+QBLKSZ-1)/QBLKSZ, (imH+QBLKSZ-1)/QBLKSZ);
	dim3 hgriddim((Rmax+QBLKSZ-1)/QBLKSZ, (angles+QBLKSZ-1)/QBLKSZ);
	CUALLOC(ima_d, sz*sizeof(float));
	CUMOV2DEV(ima_d, ima, sz*sizeof(float));
	sz = Rmax * angles;
	CUALLOC(hough_d, sz*sizeof(float));
	fill_zeros<<<hgriddim, blkdim>>>(hough_d, sz, Rmax, angles); // обнуляем
	cudaThreadSynchronize();
	//~ CUMOV2DEV(hough_d, hough, sz*sizeof(float));
	fill_lin_hough_array<<<griddim, blkdim>>>(ima_d, imW, imH, Sin_d, Cos_d,
									Rmax, angles, treshold, hough_d);
	cudaThreadSynchronize();
	normalize_vec<<<1, lblksz>>>(hough_d, sz);
	cudaThreadSynchronize();
	CUMOV2HOST(hough, hough_d, sz*sizeof(float));
free_all:
	CUFREE(hough_d);
	CUFREE(ima_d);
	return ret;
	#undef RETMACRO
	#define RETMACRO return
}