Емельянов Эдуард Владимирович (eddy_em) wrote,
Емельянов Эдуард Владимирович
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 
</pre>
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 — из-за моей ошибки):
randnum
Генерирование случайных чисел.
interplin interplog
Интерполяция. Справа — в логарифмическом масштабе.
Tags: c, cuda, обработка изображений, сплайны
Subscribe

  • Опять Подорванка смыла мост

    В четверг лило настолько, что, похоже, опять на подорванке забилась стремнина бревнами, а потом внезапно это все прорвало. МЧСовсцы перетаскивали…

  • Diskworld

    Смотрели вчера с женой на кухне "Крепкого орешка" и зацепился мой взгляд за дочкину поделку: Я жене и говорю: прямо-таки просится эти белые…

  • Юпитер и Сатурн

    Слишком близко к Солнцу, чтобы что-то получилось, а также слишком разные по яркости, чтобы без обработки их можно было бы снять одновременно.

promo eddy_em august 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic
  • 4 comments