Фильтр быстрого размытия по Гауссу с ARM NEON

Большая часть проблемы здесь в том, что алгоритм слишком точен, как указал @PaulR. Обычно лучше, чтобы ваша таблица коэффициентов была не более точной, чем ваши данные. В этом случае, поскольку кажется, что вы обрабатываете uchar данные, вы должны использовать примерно 8-битную таблицу коэффициентов.

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

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

Это могло бы сработать лучше, если бы вы использовали специальную версию цикла около краев, а затем, когда вы в безопасности от этого, вы используете упрощенный внутренний цикл, который не вызывает эту reflect101() функцию.

Во-вторых (более актуально для кода прототипа), можно сложить крылья окна вместе перед применением весовой функции, потому что таблица содержит одинаковые коэффициенты с обеих сторон.

sum = src.at<uchar>(y1, x) * coeffs[3];
for(int i = -3; i < 0; i++) {
    int tmp = src.at<uchar>(y + i, x) + src.at<uchar>(y - i, x);
    sum += coeffs[i + 3] * tmp;
}

Это сэкономит вам шесть умножений на пиксель, и это шаг к некоторым другим оптимизациям, связанным с контролем условий переполнения.

Затем есть пара других проблем, связанных с системой памяти.

Двухпроходный подход в принципе хорош, поскольку избавляет вас от необходимости выполнять много повторных вычислений. К сожалению, он может вытолкнуть полезные данные из кеша L1, что может значительно замедлить работу. Это также означает, что, когда вы записываете результат в память, вы вычисляете промежуточную сумму квантованием, что может снизить точность.

Когда вы конвертируете этот код в NEON, одна из вещей, на которой вы захотите сосредоточиться, - это попытаться сохранить рабочий набор внутри файла регистров, но не отбрасывая вычисления до того, как они будут полностью использованы.

Когда люди действительно используют два прохода, обычно промежуточные данные переносятся, то есть столбец ввода становится строкой вывода.

Это связано с тем, что ЦП действительно не любит извлекать небольшие объемы данных по нескольким строкам входного изображения. Это работает намного эффективнее (из-за того, как работает кеш), если вы соберете вместе кучу горизонтальных пикселей и отфильтруете их. Если временный буфер транспонирован, то второй проход также собирает вместе группу горизонтальных точек (которые были бы вертикальными в исходной ориентации), и он снова транспонирует свой вывод, чтобы он получился правильным.

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

Если это специально для 8-битных изображений, вам действительно не нужны коэффициенты с плавающей запятой, особенно двойная точность. Также вы не хотите использовать числа с плавающей запятой для x1, y1. Вы должны просто использовать целые числа для координат, и вы можете использовать фиксированную точку (то есть целое число) для коэффициентов, чтобы сохранить всю арифметику фильтра в целочисленной области, например

Mat myGaussianBlur(Mat src){
    Mat dst(src.rows, src.cols, CV_8UC1);
    Mat temp(src.rows, src.cols, CV_16UC1); // <<<
    int sum, x1, y1;  // <<<

    // coefficients of 1D gaussian kernel with sigma = 2
    double coeffs[] = {0.06475879783, 0.1209853623, 0.1760326634, 0.1994711402, 0.1760326634, 0.1209853623, 0.06475879783};
    int coeffs_i[7] = { 0 }; // <<<
    //Normalize coeffs
    float coeffs_sum = 0.9230247873f;
    for (int i = 0; i < 7; i++){
        coeffs_i[i] = (int)(coeffs[i] / coeffs_sum * 256); // <<<
    }

    // filter vertically
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0; // <<<
            for(int i = -3; i <= 3; i++){
                y1 = reflect101(src.rows, y - i);
                sum += coeffs_i[i + 3]*src.at<uchar>(y1, x); // <<<
            }
            temp.at<uchar>(y,x) = sum;
        }
    }

    // filter horizontally
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0; // <<<
            for(int i = -3; i <= 3; i++){
                x1 = reflect101(src.rows, x - i);
                sum += coeffs_i[i + 3]*temp.at<uchar>(y, x1); // <<<
            }
            dst.at<uchar>(y,x) = sum / (256 * 256); // <<<
        }
    }

    return dst;
}

Это код после реализации всех предложений @Paul R и @ sh1, кратко изложенный следующим образом:

1) используйте только целочисленную арифметику (с точностью по вкусу)

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

3) применять только горизонтальные фильтры, чтобы воспользоваться хранилищем по строкам матриц

4) отделите циклы по краям от тех, которые находятся внутри изображения, чтобы не вызывать ненужных вызовов функций отражения. Я полностью убрал функции отражения, включив их внутри петель по краям.

5) Кроме того, в качестве личного наблюдения, чтобы улучшить округление без вызова (медленной) функции «round» или «cvRound», я добавил как к временным, так и к конечным результатам пикселей 0,5f (= 32768 с точностью до целых чисел), чтобы уменьшить ошибка / разница по сравнению с OpenCV.

Теперь производительность намного лучше, примерно в 15-6 раз ниже, чем у OpenCV.

Однако результирующая матрица не полностью идентична матрице, полученной с помощью размытия по Гауссу в OpenCV. Это связано не с арифметической длиной (достаточной), так как с устранением ошибки остается. Обратите внимание, что это минимальная разница между 0 и 2 (по абсолютной величине) интенсивности пикселей между матрицами, полученными из двух версий. Коэффициенты такие же, как и в OpenCV, полученные с помощью getGaussianKernel с тем же размером и сигмой.

Mat myGaussianBlur(Mat src){

Mat dst(src.rows, src.cols, CV_8UC1);
Mat temp(src.rows, src.cols, CV_8UC1);
int sum;
int x1;

double coeffs[] = {0.070159, 0.131075, 0.190713, 0.216106, 0.190713, 0.131075, 0.070159};
int coeffs_i[7] = { 0 };
for (int i = 0; i < 7; i++){
        coeffs_i[i] = (int)(coeffs[i] * 65536); //65536
}

// filter horizontally - inside the image
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = 3; x < (src.cols - 3); x++){
        sum = ptr[x] * coeffs_i[3];
        for(int i = -3; i < 0; i++){
            int tmp = ptr[x+i] + ptr[x-i];
            sum += coeffs_i[i + 3]*tmp;
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
// filter horizontally - edges - needs reflect
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = 0; x <= 2; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 < 0){
                x1 = -x1;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = (src.cols - 3); x < src.cols; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 >= src.cols){
                x1 = 2*src.cols - x1 - 2;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}

// transpose to apply again horizontal filter - better cache data locality
transpose(temp, temp);

// filter horizontally - inside the image
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = 3; x < (src.cols - 3); x++){
        sum = ptr[x] * coeffs_i[3];
        for(int i = -3; i < 0; i++){
            int tmp = ptr[x+i] + ptr[x-i];
            sum += coeffs_i[i + 3]*tmp;
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
// filter horizontally - edges - needs reflect
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = 0; x <= 2; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 < 0){
                x1 = -x1;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = (src.cols - 3); x < src.cols; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 >= src.cols){
                x1 = 2*src.cols - x1 - 2;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}

transpose(dst, dst);

return dst;
}

Согласно документу Google, на устройстве Android использование float / double в два раза медленнее, чем использование int / uchar.

Вы можете найти некоторые решения для ускорения кода C ++ в этих документах Android. https://developer.android.com/training/articles/perf-tips