Реализации алгоритмов/Билинейная интерполяция

Билинейная интерполяция — в вычислительной математике расширение линейной интерполяции для функций двух переменных. Ключевая идея заключается в том, чтобы провести обычную линейную интерполяцию сначала в одном направлении, затем в перпендикулярном. Формула билинейной интерполяции интерполирует значения функции в произвольном прямоугольнике по четырем её значениям в вершинах прямоугольника и экстраполирует функцию на всю остальную плоскость.

Ниже приведен пример программы билинейной интерполяции изображения, написанный на C99 (на C89 компилироваться не будет!)

Входные параметры:

  a       - указатель на массив пикселей изображения, которое необходимо увеличить (уменьшить)
            Нумерация элементов [0..old_h-1, 0..old_w-1]
  oldw   - старая ширина изображения
  oldh   - старая высота изображения

Выходные параметры:

   b      - указатель на массив пикселей ресемплированного изображения
            Нумерация элементов [0..new_h-1, 0..new_w-1]
   neww  - новая ширина изображения
   newh  - новая высота изображения
#include <stdio.h>
#include <math.h>
#include <sys/types.h>

void resample(int oldw, int oldh, int neww, int newh, u_int a[oldh][oldw], u_int b[newh][neww])
{
	int i, j;
	int h, w;
	float t;
	float u;
	float tmp;
	float d1, d2, d3, d4;
	u_int p1, p2, p3, p4;	/* Окрестные пикселы */

	u_char red, green, blue;

	for (j = 0; j < newh; j++) {
		tmp = (float) (j) / (float) (newh - 1) * (oldh - 1);
		h = (int) floor(tmp);
		if (h < 0) {
			h = 0;
		} else {
			if (h >= oldh - 1) {
				h = oldh - 2;
			}
		}
		u = tmp - h;

		for (i = 0; i < neww; i++) {

			tmp = (float) (i) / (float) (neww - 1) * (oldw - 1);
			w = (int) floor(tmp);
			if (w < 0) {
				w = 0;
			} else {
				if (w >= oldw - 1) {
					w = oldw - 2;
				}
			}
			t = tmp - w;

			/* Коэффициенты */
			d1 = (1 - t) * (1 - u);
			d2 = t * (1 - u);
			d3 = t * u;
			d4 = (1 - t) * u;

			/* Окрестные пиксели: a[i][j] */
			p1 = a[h][w];
			p2 = a[h][w + 1];
			p3 = a[h + 1][w + 1];
			p4 = a[h + 1][w];

			/* Компоненты */
			blue = (u_char) p1 *d1 + (u_char) p2 *d2 + (u_char) p3 *d3 + (u_char) p4 *d4;
			green = (u_char) (p1 >> 8) * d1 + (u_char) (p2 >> 8) * d2 + (u_char) (p3 >> 8) * d3 + (u_char) (p4 >> 8) * d4;
			red = (u_char) (p1 >> 16) * d1 + (u_char) (p2 >> 16) * d2 + (u_char) (p3 >> 16) * d3 + (u_char) (p4 >> 16) * d4;

			/* Новый пиксел из R G B  */
			b[j][i] = ((u_int32_t) red << 16) | ((u_int32_t) green << 8) | (blue);
		}
	}
}