Оптимизация вычисления интегральной матрицы

У меня есть одноканальная матрица размерности width x height. Мне необходимо вычислить на ее основе интегральную матрицу, где каждый ее элемент на позиции [x, y] будет равен сумме всех элементов изначальной матрицы в прямоугольнике [0, 0, x, y]. Простейший алгоритм. Как его ускорить? Очевидный вариант: разбить вычисление на два этапа — сначала вычисляем суммы по горизонтали и записываем в выходную матрицу, а потом вычисляем суммы по вертикали на основе предыдущего шага. Каждую горизонталь и вертикаль обрабатываем в отдельном потоке и в результате все должно происходить очень быстро. Однако когда я заимплементил это на cuda (два ядра: одно для горизонтали, второе для вертикали) и сравнил результат с cv::integral, то был неприятно удивлен - cv::integral работало в 4-5 раз быстрее, чем вычисления на cuda (я засекал именно вычисления, копирование с cpu на gpu и обратно не производилось)! Значит, в opencv используется какой-то гораздо более продвинутый алгоритм. Какие еще есть варианты оптимизации вычисления интегральной матрицы?