Короче говоря: я разрабатываю приложение для обработки изображений с интенсивным вычислением на С++. Он должен рассчитать множество вариантов искажений изображения на небольших блоках пикселей, извлеченных из больших изображений. Программа работает не так быстро, как хотелось бы. Профилирование (OProfile) показало, что функция деформирования/интерполяции потребляет более 70% процессорного времени, поэтому казалось, что нужно попробовать и оптимизировать ее.
Я использовал библиотеку обработки изображений OpenCV для задачи до сих пор:
// some parameters for the image warps (position, stretch, skew)
struct WarpParams;
void Image::get(const WarpParams ¶ms)
{
// fills matrices mapX_ and mapY_ with x and y coordinates of points to be
// inteprolated.
updateCoordMaps(params);
// perform interpolation to obtain pixels at point locations
// P(mapX_[i], mapY_[i]) based on the original data and put the
// result in pixels_. Use bicubic inteprolation.
cv::remap(image_->data(), pixels_, mapX_, mapY_, CV_INTER_CUBIC);
}
Я написал свою собственную функцию интерполяции и поместил ее в тестовую проводку, чтобы обеспечить правильность, пока я экспериментирую и сравниваю ее со старым.
Моя функция работала очень медленно, что и следовало ожидать. Как правило, идея заключается в следующем:
- Итерации над картами X_, mapY_ координатными картами, извлечение (вещественные) координаты следующего пикселя для интерполяции;
- Получить 4x4 соседний пиксельных значений (целочисленные координаты) из исходного изображения, окружающего интерполированный пиксель;
- Вычислить коэффициенты ядра свертки для каждого из этих 16 пикселей;
- Вычислить значение интерполированного пикселя как линейную комбинацию из 16 значений пикселей и коэффициентов ядра.
Старая функция приурочена к 25us на моем Wolfdale Core2 Duo. Новый взял 587us (!). Я с готовностью надел волшебную шляпу и начал взламывать код. Мне удалось удалить все ветки, опустить некоторые дублирующие вычисления и преобразовать 3 вложенных цикла в одну по координатным картам. Вот что я придумал:
void Image::getConvolve(const WarpParams ¶ms)
{
__declspec(align(16)) static float kernelX[4], kernelY[4];
// grab pointers to coordinate map matrices and the original image
const float
*const mapX = mapX_.ptr<float>(),
*const mapY = mapY_.ptr<float>(),
*const img = image_->data().ptr<float>();
// grab pointer to the output image
float *const subset = pixels_.ptr<float>(),
x, y, xint, yint;
const ptrdiff_t imgw = image_->width();
ptrdiff_t imgoffs;
__m128 v_px, v_kernX, v_kernY, v_val;
// iterate over the coordinate matrices as linear buffers
for (size_t idx = 0; idx < pxCount; ++idx)
{
// retrieve coordinates of next pixel from precalculated maps,
// break up each into fractional and integer part
x = modf(mapX[idx], &xint);
y = modf(mapY[idx], &yint);
// obtain offset of the top left pixel from the required 4x4
// neighborhood of the current pixel in the image
// buffer (sadly, the position will be unaligned)
imgoffs = (((ptrdiff_t)yint - 1) * imgw) + (ptrdiff_t)xint - 1;
// calculate all 4 convolution kernel values for every row and
// every column
tap4Kernel(x, kernelX);
tap4Kernel(y, kernelY);
// load the kernel values for the columns, these don't change
v_kernX = _mm_load_ps(kernelX);
// process a row of the 4x4 neighborhood
// get set of convolution kernel values for the current row
v_kernY = _mm_set_ps1(kernelY[0]);
v_px = _mm_loadu_ps(img + imgoffs); // load the pixel values
// calculate the linear combination of the pixels with kernelX
v_px = _mm_mul_ps(v_px, v_kernX);
v_px = _mm_mul_ps(v_px, v_kernY); // and kernel Y
v_val = v_px; // add result to the final value
imgoffs += imgw;
// offset points now to next row of the 4x4 neighborhood
v_kernY = _mm_set_ps1(kernelY[1]);
v_px = _mm_loadu_ps(img + imgoffs);
v_px = _mm_mul_ps(v_px, v_kernX);
v_px = _mm_mul_ps(v_px, v_kernY);
v_val = _mm_add_ps(v_val, v_px);
imgoffs += imgw;
/*... same for kernelY[2] and kernelY[3]... */
// store resulting interpolated pixel value in the subset's
// pixel matrix
subset[idx] = horizSum(v_val);
}
}
// Calculate all 4 values of the 4-tap convolution kernel for 4 neighbors
// of a pixel and store them in an array. Ugly but fast.
// The "arg" parameter is the fractional part of a pixel coordinate, i.e.
// a number in the range <0,1)
void Image::tap4Kernel(const float arg, float *out)
{
// chaining intrinsics was slower, so this is done in separate steps
// load the argument into 4 cells of a XMM register
__m128
v_arg = _mm_set_ps1(arg),
v_coeff = _mm_set_ps(2.0f, 1.0f, 0.0f, -1.0f);
// subtract vector of [-1, 0, 1, 2] to obtain coorinates of 4 neighbors
// for kernel calculation
v_arg = _mm_sub_ps(v_arg, v_coeff);
// clear sign bits, this is equivalent to fabs() on all 4
v_coeff = _mm_set_ps1(-0.f);
v_arg = _mm_andnot_ps(v_coeff, v_arg);
// calculate values of abs(argument)^3 and ^2
__m128
v_arg2 = _mm_mul_ps(v_arg, v_arg),
v_arg3 = _mm_mul_ps(v_arg2, v_arg),
v_val, v_temp;
// calculate the 4 kernel values as
// arg^3 * A + arg^2 * B + arg * C + D, using
// (A,B,C,D) = (-0.5, 2.5, -4, 2) for the outside pixels and
// (1.5, -2.5, 0, 1) for inside
v_coeff = _mm_set_ps(-0.5f, 1.5f, 1.5f, -0.5f);
v_val = _mm_mul_ps(v_coeff, v_arg3);
v_coeff = _mm_set_ps(2.5f, -2.5f, -2.5f, 2.5f);
v_temp = _mm_mul_ps(v_coeff, v_arg2);
v_val = _mm_add_ps(v_val, v_temp);
v_coeff = _mm_set_ps(-4.0f, 0.0f, 0.0f, -4.0f),
v_temp = _mm_mul_ps(v_coeff, v_arg);
v_val = _mm_add_ps(v_val, v_temp);
v_coeff = _mm_set_ps(2.0f, 1.0f, 1.0f, 2.0f);
v_val = _mm_add_ps(v_val, v_coeff);
_mm_store_ps(out, v_val);
}
Мне было приятно, что время работы на этом было ниже 40us, даже до введения SSE в основной цикл, который я сохранил для последнего. Я ожидал по крайней мере 3-кратного ускорения, но он работал только чуть быстрее на 36us, медленнее, чем старый get(), который я пытался улучшить. Хуже всего то, что, когда я изменил бенчмарк, чтобы делать больше прогонов, у старой функции было такое же среднее время выполнения, в то время как моя растянута до более 127us, а это означает, что для некоторых крайних значений параметра warp требуется еще больше времени (имеет смысл, потому что больше warps означает, что мне нужно достичь широкодисперсных значений пикселей из исходного изображения, чтобы вычислить результат).
Я понял, что причиной должны быть неравнозначные нагрузки, но это не может помочь (мне нужно достичь непредсказуемых значений пикселей). Я не мог видеть ничего больше, что мог сделать в отделе оптимизации, поэтому решил посмотреть на функцию cv:: remap(), чтобы увидеть, как они это делают. Представьте мое удивление в том, что он содержит беспорядок вложенных петель и множество ветвей. Они также делают много аргументов, которые мне не нужно беспокоиться. Насколько я могу судить (без комментариев в коде), SSE (с невыложенными нагрузками также!) Используется только для извлечения значений из карт координат и округления их в целые числа, затем вызывается функция, которая делает фактическую интерполяцию с регулярной флоат-арифметикой.
Мой вопрос: почему я так неудачно потерпел неудачу (почему мой код настолько медленный, и почему их быстрее, хотя это похоже на беспорядок) и что я могу сделать, чтобы улучшить свой код?
Я не вставляю код OpenCV здесь, потому что это уже слишком долго, вы можете проверить его на pastebin.
Я тестировал и скомпилировал свой код в режиме выпуска в VС++ 2010. Используемый OpenCV представлял собой прекомпилированный бинарный пакет v2.3.1.
EDIT: Значения пикселей - это float в диапазоне 0..1. Профилирование показало, что функция tap4Kernel() не имеет значения, большинство времени проводится внутри getConvolve().
EDIT2: Я вставил разбор сгенерированного кода в pastebin. Это скомпилировано на старом процессоре Banias Celeron (имеет SSE2), но выглядит более или менее одинаковым.
EDIT3: После чтения Что каждый программист должен знать о памяти Я понял, что неправильно предполагал, что функция OpenCV реализует более или менее одинаковые как я это делал, чего не должно быть. Для каждого пикселя, который я интерполирую, мне нужно получить его 4x4 окрестности, пиксели которых не помещаются внутри буфера изображения. Я злоупотребляю кэшами процессора, и OpenCV, вероятно, этого не делает. Профилирование VTune, похоже, соглашается с тем, что моя функция имеет 5 800 000 обращений к памяти, в то время как OpenCV выполняет только 400 000. Их функция беспорядок и, вероятно, может быть дополнительно оптимизирована, но ей все же удается преодолеть меня, возможно, из-за более разумного подхода к использованию памяти и кеша.
UPDATE:Мне удалось улучшить способ загрузки значений пикселей в регистры XMM. Я выделяю буфер в объекте, который содержит ячейки из 16 элементов для каждого пикселя изображения. При загрузке изображения я заполняю этот буфер ячейки предварительно упорядоченными последовательностями 4x4-окрестностей для каждого пикселя. Не очень экономично (изображение занимает 16 раз), но в этом случае нагрузки всегда выравниваются (не больше _mm_loadu_ps()), и я избегаю необходимости делать разбросанные чтения пикселей из буфера изображения, поскольку требуемые пиксели сохраняются последовательно. К моему удивлению, никакого улучшения вообще не было. Я слышал, что неуравновешенные нагрузки могут быть на 10 раз медленнее, но, очевидно, это не проблема. Но, комментируя части кода, я узнал, что вызовы modf() отвечают за 75% времени выполнения! Я сосредоточусь на устранении этих проблем и отправке ответов.
