eddy_em

Только я пожаловался, что давненько не писал, как — вуаля. Решил поделиться. Итак, проверил я алгоритм, про который в предыдущей заметке писал (по статье Ruijters, Romeny & Suetens «Efficient GPU-Based Texture Interpolation using Uniform B-Splines»). Получается действительно черт знает что:

idata-bad — Плохой алгоритм. Сверху — исходные данные, Снизу — интерполяция.

odata-bad — Плохой алгоритм. Сверху — исходные данные, Снизу — интерполяция.

А вот — сами тестовые данные. Было:

0.669 2.691 3.733 4.545 4.625 
1.856 1.570 2.087 4.991 5.501 
0.629 2.275 3.669 3.256 5.650 
1.420 2.747 3.031 3.467 5.347 
1.894 2.683 3.470 3.712 5.486

Стало:

0.669 0.669 0.669 1.155 1.443 1.760 2.067 2.327 2.976 3.130 3.294 3.448 3.880 3.984 4.107 4.234 4.350 4.562 4.573 4.585 
0.669 0.669 0.669 1.155 1.443 1.760 2.067 2.327 2.976 3.130 3.294 3.448 3.880 3.984 4.107 4.234 4.350 4.562 4.573 4.585 
0.669 0.669 0.669 1.155 1.443 1.760 2.067 2.327 2.976 3.130 3.294 3.448 3.880 3.984 4.107 4.234 4.350 4.562 4.573 4.585 
0.955 0.955 0.955 1.307 1.516 1.746 1.969 2.157 2.672 2.808 2.952 3.087 3.575 3.744 3.943 4.149 4.336 4.691 4.716 4.744 
1.124 1.124 1.124 1.397 1.559 1.738 1.910 2.057 2.492 2.616 2.749 2.873 3.394 3.601 3.846 4.098 4.328 4.767 4.800 4.838 
1.310 1.310 1.310 1.496 1.607 1.729 1.846 1.946 2.294 2.406 2.526 2.638 3.195 3.444 3.739 4.043 4.319 4.851 4.893 4.942 
1.490 1.490 1.490 1.592 1.653 1.720 1.784 1.839 2.102 2.202 2.309 2.410 3.002 3.292 3.635 3.989 4.311 4.932 4.984 5.042 
1.642 1.642 1.642 1.673 1.692 1.712 1.732 1.748 1.939 2.030 2.126 2.216 2.839 3.164 3.548 3.943 4.303 5.001 5.060 5.127 
1.520 1.520 1.520 1.579 1.613 1.651 1.688 1.719 1.970 2.082 2.201 2.313 2.880 3.136 3.439 3.752 4.036 4.731 4.870 5.029 
1.339 1.339 1.339 1.466 1.541 1.624 1.704 1.772 2.110 2.241 2.380 2.511 3.025 3.219 3.448 3.683 3.898 4.532 4.710 4.912 
1.146 1.146 1.146 1.346 1.465 1.595 1.721 1.828 2.258 2.410 2.571 2.722 3.180 3.306 3.456 3.610 3.751 4.322 4.540 4.787 
0.965 0.965 0.965 1.233 1.393 1.567 1.737 1.881 2.398 2.569 2.750 2.920 3.325 3.389 3.464 3.541 3.612 4.124 4.379 4.670 
0.771 0.771 0.771 1.153 1.380 1.628 1.870 2.074 2.687 2.863 3.051 3.227 3.507 3.474 3.434 3.394 3.357 3.775 4.089 4.445 
0.873 0.873 0.873 1.245 1.466 1.708 1.943 2.142 2.709 2.864 3.029 3.184 3.445 3.426 3.403 3.379 3.357 3.789 4.093 4.439 
0.993 0.993 0.993 1.353 1.567 1.802 2.030 2.222 2.734 2.865 3.003 3.134 3.371 3.368 3.365 3.362 3.359 3.804 4.098 4.432 
1.117 1.117 1.117 1.465 1.672 1.899 2.119 2.305 2.760 2.865 2.977 3.081 3.295 3.310 3.326 3.344 3.360 3.821 4.103 4.425 
1.230 1.230 1.230 1.567 1.767 1.987 2.200 2.381 2.784 2.866 2.952 3.034 3.226 3.256 3.291 3.328 3.361 3.835 4.108 4.418 
1.519 1.519 1.519 1.811 1.984 2.174 2.359 2.515 2.840 2.898 2.959 3.016 3.194 3.245 3.305 3.367 3.423 3.907 4.160 4.447 
1.583 1.583 1.583 1.858 2.021 2.199 2.373 2.519 2.850 2.918 2.990 3.057 3.249 3.296 3.352 3.410 3.463 3.937 4.188 4.473 
1.657 1.657 1.657 1.911 2.062 2.228 2.388 2.525 2.862 2.941 3.025 3.104 3.311 3.355 3.406 3.460 3.508 3.971 4.220 4.503

Ясно, что такое уродство ну совсем не годится для интерполяции данных: более наглядно тогда уж просто «родными» средствами GPU сделать интерполяцию текстуры. Зашел я в википедию и почитал все-таки статью о бикубической интерполяции. Расписал на бумажечке простейшее матричное выражение и получил вот такой код:


inline __device__ float p_ta(float t, float t2, float t3,
						float a,float b,float c,float d){
	return (2*b + t*(-a+c) + t2*(2*a-5*b+4*c-d) + t3*(-a+3*b-3*c+d)) / 2.f;
}

__device__ float interpolate_bicubic(float x, float y){
	#define TEX(x,y) tex2D(data_texture, x, y)
	#define TX(Y) p_ta(pt.x, pt2.x, pt3.x, TEX(x0-1,y0+(Y)), TEX(x0,y0+(Y)), TEX(x0+1,y0+(Y)), TEX(x0+2,y0+(Y)))
	float2 crd = make_float2(x,y);
	float2 zeropt = floor(crd);
	float2 pt = crd - zeropt;
	float2 pt2 = pt*pt;
	float2 pt3 = pt*pt2;
	float x0 = zeropt.x, y0 = zeropt.y;
	//return x;
	return p_ta(
		pt.y, pt2.y, pt3.y,
		TX(-1), TX(0), TX(1), TX(2)
	);
	#undef TEX
	#undef TX
}
__global__ void bicubic_interp(float *out,
								float fracX, float fracY,
								unsigned int oW, unsigned int oH){
	int X, Y;  // pixel coordinates on output
	float x,y; // coordinates on output in value of input
	X = threadIdx.x + blockDim.x * blockIdx.x;
	Y = threadIdx.y + blockDim.y * blockIdx.y;
	if(X >= oW || Y >= oH) return;
	x = (float)X * fracX;
	y = (float)Y * fracY;
	//out[Y*oW+X] = interpolate_bspline(x, y);
	out[Y*oW+X] = interpolate_bicubic(x, y);
}

EXTERN int CUbicubic_interp(float *out, float *in,
								size_t oH, size_t oW, size_t iH, size_t iW){
	FNAME();
	float *in_dev = NULL, *out_dev = NULL;
	size_t iWp = iW*sizeof(float), oSz = oH*oW*sizeof(float), pitch;
	ret = 1;
	dim3 blkdim(QBLKSZ, QBLKSZ);
	dim3 griddim((oW+QBLKSZ-1)/QBLKSZ, (oH+QBLKSZ-1)/QBLKSZ);
	CUALLOCPITCH(in_dev, &pitch, iWp, iH);
	CUMOV2DEVPITCH(in_dev, pitch, in, iWp, iH);
	CUALLOC(out_dev, oSz);

	CUerr = cudaBindTexture2D(NULL, data_texture, in_dev, iW, iH, pitch);
	data_texture.addressMode[0] = cudaAddressModeClamp;
	data_texture.addressMode[1] = cudaAddressModeClamp;
	data_texture.normalized = false;  // access with absolute texture coordinates
	data_texture.filterMode = cudaFilterModePoint; // cudaFilterModeLinear
	if(CUERROR("CUDA: can't bind texture")){FREERETMACRO;}

	bicubic_interp<<<griddim, blkdim>>>(out_dev, (float)(iW-1)/(oW-1), (float)(iH-1)/(oH-1), oW, oH);
	cudaThreadSynchronize();
	CUMOV2HOST(out, out_dev, oSz);
free_all:
	CUFREE(in_dev);
	CUFREE(out_dev);
	return ret;
}

Запись самой __device__-функции получилась даже короче, чем у «упрощенного B-сплайна». И результаты значительно приличнее:

idata — Хороший алгоритм. Слева — исходные данные, справа — интерполяция.

odata — Хороший алгоритм. Слева — исходные данные, справа — интерполяция.

А вот и сами тестовые данные. Было:

0.760 1.069 3.221 4.291 4.968 
1.716 2.277 3.396 3.561 4.659 
1.172 2.280 3.989 3.936 5.978 
0.158 2.035 2.001 4.650 4.540 
1.705 1.778 3.728 4.701 4.435

Стало:

0.760 0.773 0.774 0.807 0.914 1.140 1.532 2.032 2.559 3.031 3.382 3.653 3.874 4.065 4.245 4.427 4.598 4.751 4.877 4.968 
0.925 0.946 0.960 1.005 1.117 1.332 1.695 2.155 2.638 3.068 3.382 3.616 3.802 3.963 4.123 4.300 4.480 4.650 4.793 4.895 
1.172 1.204 1.238 1.300 1.416 1.614 1.934 2.334 2.750 3.116 3.378 3.558 3.692 3.809 3.940 4.108 4.301 4.493 4.659 4.776 
1.433 1.479 1.535 1.617 1.740 1.922 2.197 2.534 2.879 3.178 3.382 3.503 3.579 3.648 3.747 3.907 4.115 4.335 4.531 4.666 
1.639 1.700 1.779 1.883 2.018 2.189 2.431 2.719 3.008 3.252 3.407 3.475 3.497 3.521 3.592 3.750 3.980 4.233 4.463 4.620 
1.723 1.799 1.900 2.027 2.178 2.354 2.586 2.857 3.123 3.342 3.471 3.501 3.483 3.473 3.525 3.690 3.949 4.243 4.511 4.693 
1.692 1.781 1.902 2.052 2.226 2.423 2.684 2.984 3.276 3.509 3.631 3.626 3.560 3.508 3.542 3.728 4.040 4.402 4.734 4.960 
1.584 1.687 1.828 2.001 2.204 2.435 2.743 3.101 3.446 3.715 3.843 3.807 3.692 3.596 3.613 3.828 4.207 4.651 5.062 5.340 
1.429 1.548 1.711 1.911 2.145 2.410 2.766 3.177 3.573 3.881 4.022 3.971 3.828 3.706 3.718 3.962 4.399 4.913 5.388 5.710 
1.256 1.395 1.585 1.817 2.082 2.374 2.752 3.185 3.601 3.927 4.086 4.048 3.917 3.808 3.835 4.101 4.565 5.108 5.608 5.947 
1.072 1.238 1.469 1.744 2.041 2.342 2.696 3.085 3.458 3.761 3.937 3.959 3.907 3.875 3.953 4.225 4.664 5.165 5.624 5.936 
0.799 1.009 1.309 1.653 1.994 2.291 2.566 2.834 3.088 3.320 3.519 3.665 3.788 3.922 4.103 4.372 4.728 5.108 5.448 5.682 
0.498 0.757 1.133 1.553 1.941 2.225 2.397 2.509 2.610 2.751 2.977 3.279 3.620 3.962 4.269 4.527 4.769 4.987 5.166 5.297 
0.254 0.551 0.987 1.465 1.886 2.154 2.231 2.202 2.167 2.228 2.480 2.928 3.472 4.006 4.426 4.667 4.793 4.848 4.868 4.894 
0.152 0.460 0.915 1.410 1.837 2.089 2.112 2.007 1.901 1.923 2.198 2.740 3.411 4.064 4.550 4.769 4.806 4.738 4.640 4.588 
0.275 0.552 0.961 1.406 1.794 2.033 2.069 1.994 1.924 1.976 2.267 2.820 3.497 4.148 4.622 4.814 4.804 4.681 4.534 4.448 
0.606 0.812 1.112 1.446 1.753 1.974 2.077 2.120 2.175 2.314 2.617 3.116 3.705 4.260 4.659 4.811 4.779 4.645 4.491 4.399 
1.031 1.149 1.313 1.508 1.715 1.919 2.114 2.315 2.535 2.789 3.100 3.514 3.967 4.379 4.672 4.782 4.742 4.621 4.486 4.404 
1.436 1.470 1.507 1.570 1.684 1.871 2.157 2.510 2.889 3.251 3.569 3.897 4.214 4.485 4.675 4.746 4.706 4.605 4.495 4.427 
1.705 1.683 1.634 1.610 1.661 1.838 2.182 2.636 3.120 3.554 3.877 4.150 4.379 4.558 4.680 4.726 4.683 4.592 4.496 4.435

UPD:
Сделал реализацию на CPU. Долго ломал голову, как бы оптимизировать, но в конце-концов плюнул и сделал тупо:


inline float p_ta(float t, float t2, float t3,
						float a,float b,float c,float d){
	return (2*b + t*(-a+c) + t2*(2*a-5*b+4*c-d) + t3*(-a+3*b-3*c+d)) / 2.f;
}

int CPUbicubic_interp(float *out, float *in,
								size_t oH, size_t oW, size_t iH, size_t iW){
	FNAME();
	float fracX = (float)(iW-1)/(oW-1), fracY = (float)(iH-1)/(oH-1);
	size_t X, Y, ym1,y1,y2, xm1,x1,x2;  // pixel coordinates on output
	size_t Ym = iH - 1, Xm = iW - 1, Pcur;
	float x,y; // coordinates on output in value of input
	OMP_FOR
	for(Y = 0; Y < oH; Y++){
		// we can't do "y+=fracY" because of possible cumulative error
		y = (float)Y * fracY;
		int y0 = floor(y);
		float pty = y - (float)y0, pty2 = pty*pty, pty3 = pty*pty2;
		ym1 = y0-1; y1 = y0 + 1; y2 = y0 + 2;
		if(y0 == 0) ym1 = 0;
		if(y1 > Ym) y1 = Ym;
		if(y2 > Ym) y2 = Ym;
		for(X = 0, Pcur = Y * oW; X < oW; X++, Pcur++){
			// we can't do "x+=fracX" because of possible cumulative error
			x = (float)X * fracX;
			int x0 = floor(x);
			float ptx = x - (float)x0, ptx2 = ptx*ptx, ptx3 = ptx*ptx2;
			xm1 = x0-1; x1 = x0 + 1; x2 = x0 + 2;
			if(x0 == 0) xm1 = 0;
			if(x1 > Xm) x1 = Xm;
			if(x2 > Xm) x2 = Xm;
			#define TEX(x,y) (in[iW*y + x])
			#define TX(y) p_ta(ptx, ptx2, ptx3, TEX(xm1,y), \
						TEX(x0,y), TEX(x1,y), TEX(x2,y))
			out[Pcur] = p_ta(
				pty, pty2, pty3,
				TX(ym1), TX(y0), TX(y1), TX(y2));
			#undef TX
			#undef TEX
		}
	}
	return 1;
}

Недостаток оптимизации хорошо отразился на времени выполнения. Тесты на время выполнения (завал на графике ГСЧ на GPU — из-за моей ошибки):