Двумерная интерполяция кубическими сплайнами

[eddy_em]

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

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

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

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-сплайна». И результаты значительно приличнее:

 

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

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

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 — из-за моей ошибки):

Генерирование случайных чисел.

 

Интерполяция. Справа — в логарифмическом масштабе.

Comments

[eddy-em.livejournal.com]

Графики наспех построил. В принципе, особой пользы от них нет: просто отражают то, что в куде ГСЧ работает с постоянной скоростью, независимо от объема данных, а вот интерполяция даже на не очень крутом GPU выполняется на порядок быстрее, чем на GPU, что важно для сложных вычислений.

Как закончу — выложу, как обычно, код на сосфорж. У меня разработка затягивается раза в три еще и из-за того, что я решил дублировать код на CUDA кодом с использованием OpenMP. Без текстур сложнее параллелить…

А насчет B-сплайнов хорошо в википедии (английской) написано. Зачем дублировать?

Я, кстати, еще не до конца придумал алгоритм трассировки лучей: с одной стороны вычисления хочется ускорить, а с другой — лень много оптимизировать. А то ведь можно до того дооптимизироваться, что хрен что работать будет!