opencv - 将 OpenCV 的 Mat 容器与用于矩阵乘法的 blas 接口(interface)

Question

我正在处理超高清 (2160 x 3840) 图像。我所做的其中一项处理包括在 X 轴和 Y 轴上处理 Sobel 滤波，然后我必须将每个输出矩阵乘以它的转置，然后我将梯度图像处理为梯度总和的平方根。

所以:S = sqrt( S_x * S_x^t + S_y * S_y^t)。

由于图像的尺寸，OpenCV 在不使用多线程的情况下最多需要 20 秒来处理该图像，而使用多线程则需要 10 秒。

我知道 OpenCV 会调用 OpenCL 以加快过滤操作，因此我认为可能需要很长时间才能尝试从过滤步骤中获得性能。

对于矩阵乘法，我从 OpenCV 的 OpenCL gemm 内核实现中体验到一种不稳定性。

所以我想尝试使用 OpenBLAS insted。

我的问题是:

1.)

我编写了以下代码，但我遇到了 OpenCV 的 Mat 对象接口(interface)的一些问题:

template<class _Ty>
void mm(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    static_assert(true,"support matrix_multiply is only defined for floating precision numbers.");
}

template<>
inline void mm<float>(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    const int M = A.rows;
    const int N = B.cols;
    const int K = A.cols;

    cblas_sgemm( CblasRowMajor ,// 1
                 CblasNoTrans, // 2 TRANSA
                 CblasNoTrans, // 3 TRANSB
                 M,       // 4 M
                 N,       // 5 N
                 K,       // 6 K
                 1.,           // 7 ALPHA
                 A.ptr<float>(),//8 A
                 A.rows,        //9 LDA
                 B.ptr<float>(),//10 B
                 B.rows,        //11 LDB
                 0.,            //12 BETA
                 C.ptr<float>(),//13 C
                 C.rows);       //14 LDC

}

template<>
inline void mm<double>(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,A.rows,B.cols,A.cols,1.,A.ptr<double>(),A.rows,B.ptr<double>(),B.cols,0.,C.ptr<double>(),C.rows);
}
    void matrix_multiply(cv::InputArray _src1, cv::InputArray _src2, cv::OutputArray _dst)
    {

        CV_DbgAssert(  (_src1.isMat() || _src1.isUMat()) && (_src1.kind() == _src2.kind()) &&
                      (_src1.depth() == _src2.depth()) && (_src1.depth() == CV_32F) && (_src1.depth() == _src1.type()) &&
                      (_src1.rows() == _src2.cols())
                       );


        cv::Mat src1 = _src1.getMat();
        cv::Mat src2 = _src2.getMat();
        cv::Mat dst;

        bool cpy(false);

        if(_dst.rows() == _src1.rows() && _dst.cols() == _src2.cols() && _dst.type() == _src1.type())
            dst = _dst.getMat();
        else
        {
            dst = cv::Mat::zeros(src1.rows,src2.cols,src1.type());
            cpy = true;
        }

        if(cpy)
            dst.copyTo(_dst);
    }

我尝试按照此处指定的方式组织数据: http://www.netlib.org/lapack/explore-html/db/dc9/group__single__blas__level3.html#gafe51bacb54592ff5de056acabd83c260

没有成功。 这是我的主要问题

2.) 我在想，为了尝试加快我的实现速度，以应用此处说明的分而治之方法:

https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm

但是对于只有四个子矩阵。有没有人尝试过类似的方法或获得更好的方法来获得矩阵乘法的性能(不使用 GPU)？

提前感谢您的帮助。

Answer 1

我找到了问题 1) 的解决方案。我的第一个实现基于 BLAS 库的文档。BLAS 是用 Fortran 语言编写的，在这种语言中，索引从 1 开始，而不是像 C 或 C++ 中那样从 0 开始。另一件事是许多用 Fortran 语言编写的库按列顺序组织内存(例如 BLAS、LAPACK)，而不是大多数 C 或 C++ 库(例如 OpenCV)按行顺序组织内存。

在计算这两个属性后，我将代码修改为:

template<class _Ty>
void mm(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    static_assert(true,"The function gemm is only defined for floating precision numbers.");
}

template<>
void mm<float>(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    const int M = A.cols+1;
    const int N = B.rows;
    const int K = A.cols;

    cblas_sgemm( CblasRowMajor ,// 1
                 CblasNoTrans, // 2 TRANSA
                 CblasNoTrans, // 3 TRANSB
                 M,       // 4 M
                 N,       // 5 N
                 K,       // 6 K
                 1.,           // 7 ALPHA
                 A.ptr<float>(),//8 A
                 A.step1(),        //9 LDA
                 B.ptr<float>(),//10 B
                 B.step1(),        //11 LDB
                 0.,            //12 BETA
                 C.ptr<float>(),//13 C
                 C.step1());       //14 LDC
}

template<>
void mm<double>(cv::Mat& A,cv::Mat& B,cv::Mat& C)
{
    const int M = A.cols+1;
    const int N = B.rows;
    const int K = A.cols;

    cblas_dgemm( CblasRowMajor ,// 1
                 CblasNoTrans, // 2 TRANSA
                 CblasNoTrans, // 3 TRANSB
                 M,       // 4 M
                 N,       // 5 N
                 K,       // 6 K
                 1.,           // 7 ALPHA
                 A.ptr<double>(),//8 A
                 A.step1(),        //9 LDA
                 B.ptr<double>(),//10 B
                 B.step1(),        //11 LDB
                 0.,            //12 BETA
                 C.ptr<double>(),//13 C
                 C.step1());       //14 LDC
}

而且每件事都运作良好。在没有额外的多线程或分而治之的方法的情况下，我能够将我的代码一步的处理时间从 150 毫秒减少到 500 微秒。所以它为我解决了所有问题 :)。