Более быстрый метод множественной билинейной интерполяции?

Я пишу программу на C ++ для реконструкции 3D-объекта из набора проецируемых 2D-изображений, наиболее интенсивная часть вычислений которого заключается в увеличении и сдвиге каждого изображения с помощью билинейной интерполяции. В настоящее время у меня есть пара функций для этой задачи; «blnSetup» определяет несколько параметров вне цикла, затем «билинейный» применяет интерполяцию точка-точка внутри цикла:

(ПРИМЕЧАНИЕ: «I» — это одномерный массив, содержащий упорядоченные строки данных изображения)

//Pre-definition structure (in header)
struct blnData{
float* X;
float* Y;
int* I;
float X0;
float Y0;
float delX;
float delY;
};

//Pre-definition function (outside the FOR loop)
extern inline blnData blnSetup(float* X, float* Y, int* I)
{
blnData bln;
//Create pointers to X, Y, and I vectors
bln.X = X;
bln.Y = Y;
bln.I = I;

//Store offset and step values for X and Y
bln.X0 = X[0];
bln.delX = X[1] - X[0];
bln.Y0 = Y[0];
bln.delY = Y[1] - Y[0];

return bln;
}

//Main interpolation function (inside the FOR loop)
extern inline float bilinear(float x, float y, blnData bln)
{
float Ixy;

//Return -1 if the target point is outside the image matrix
if (x < bln.X[0] || x > bln.X[-1] || y < bln.Y[0] || y > bln.Y[-1])
Ixy = 0;
//Otherwise, apply bilinear interpolation
else
{
//Define known image width
int W = 200;

//Find nearest indices for interpolation
int i = floor((x - bln.X0) / bln.delX);
int j = floor((y - bln.Y0) / bln.delY);

//Interpolate I at (xi, yj)
Ixy = 1 / ((bln.X[i + 1] - bln.X[i])*(bln.Y[j + 1] - bln.Y[j])) *
(
bln.I[W*j + i] * (bln.X[i + 1] - x) * (bln.Y[j + 1] - y) +
bln.I[W*j + i + 1] * (x - bln.X[i]) * (bln.Y[j + 1] - y) +
bln.I[W*(j + 1) + i] * (bln.X[i + 1] - x) * (y - bln.Y[j]) +
bln.I[W*(j + 1) + i + 1] * (x - bln.X[i]) * (y - bln.Y[j])
);
}

return Ixy;
}

РЕДАКТИРОВАТЬ: вызовы функции ниже. ‘flat.imgdata’ — это std :: vector, содержащий данные входного изображения, а ‘proj.imgdata’ — это std :: vector, содержащий преобразованное изображение.

int Xs = flat.dim[0];
int Ys = flat.dim[1];

int* Iarr = flat.imgdata.data();
float II, x, y;

bln = blnSetup(X, Y, Iarr);

for (int j = 0; j < flat.imgdata.size(); j++)
{
x = 1.2*X[j % Xs];
y = 1.2*Y[j / Xs];
II = bilinear(x, y, bln);
proj.imgdata[j] = (int)II;
}

С тех пор, как я начал оптимизировать, я смог сократить время вычислений на ~ 50x (!), Переключаясь с std :: vectors на массивы C в рамках функции интерполяции, и еще на 2x или около того, очищая избыточные вычисления / typecasting / etc, но при условии, что O (n) с n является общим числом обработанных пикселей, полная реконструкция (~ 7e10 пикселей) должна все же занять 40 минут или около того — примерно на порядок больше, чем моя цель <5 минут.

Согласно профилировщику производительности Visual Studio, вызов функции интерполяции («II = bilinear (x, y, bln);»), что неудивительно, по-прежнему является основной частью моей вычислительной нагрузки. Мне не удалось найти какие-либо линейные алгебраические методы для быстрой многократной интерполяции, поэтому мой вопрос: так ли это в основном так быстро, как мой код, если не использовать больше или более быстрых процессоров для выполнения задачи? Или есть другой подход, который может ускорить процесс?

Постскриптум Я также занимаюсь кодированием на C ++ уже около месяца, поэтому не стесняйтесь указывать на любые ошибки новичка, которые я мог бы совершить.