Matlab 快速多通道积分图计算函数

所谓快速多通道积分图计算，其实就是 cumsum2D。
我写了一个比较快的版本（比 VLFeat 的快）,用 mex 编译一下就能用了。

代码

#include <string.h>
#include <mex.h>
#include <stdio.h>
#include <stdint.h>

// compute integral image
template <typename T>
void integral(const T* in, T* out, mwSize h, mwSize w, mwSize ch)
{
    mwSize iw = w + 1;
    mwSize ih = h + 1;
    out += 1;
    memset(out, 0, ih * iw * ch * sizeof(T));
    for (mwSize c = 0; c < ch; ++c)
    {
        // offset one row and one column
        out += ih;
        // from first to last column
        for (mwSize x = 0; x < w; ++x, in += h, out += ih)
        {
            T s = 0.0f;
            // from first to last column
            for (mwSize y = 0; y < h; ++y)
            {
                s += in[y];
                out[y] = out[y - ih] + s;
            }
        }
    }
}

// itg = integralChannels(channels)
void mexFunction (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
    mxClassID classId = mxGetClassID(prhs[0]);
    mwSize nDims = mxGetNumberOfDimensions(prhs[0]);
    const mwSize *dims = mxGetDimensions(prhs[0]);
    mwSize *iDims = new mwSize[nDims];
    
    if (nDims != 3 && nDims != 2)
    {
        mexErrMsgTxt("integralChannels() can only integrate single or multichannel image.");
    }
    mwSize h = dims[0];
    mwSize w = dims[1];
    mwSize ch = 1;
        
    iDims[0] = h + 1;
    iDims[1] = w + 1;
    if (nDims == 3)
    {
        ch = dims[2];
        iDims[2] = ch;
    }
    plhs[0] = mxCreateNumericArray(nDims, iDims, classId, mxREAL);
    delete [] iDims;
    switch (classId)
    {
        case mxSINGLE_CLASS:
            integral<float>((float *)mxGetData(prhs[0]), 
                    (float *)mxGetData(plhs[0]), h, w, ch);
            break;
        case mxDOUBLE_CLASS:
            integral<double>((double *)mxGetData(prhs[0]), 
                    (double *)mxGetData(plhs[0]), h, w, ch);
            break;
        case mxUINT32_CLASS:
            integral<uint32_t>((uint32_t *)mxGetData(prhs[0]), 
                    (uint32_t *)mxGetData(plhs[0]), h, w, ch);
            break;
        case mxINT32_CLASS:
            integral<int32_t>((int32_t *)mxGetData(prhs[0]), 
                    (int32_t *)mxGetData(plhs[0]), h, w, ch);
            break;
        default:
            mexErrMsgTxt("Illegal Class ID.");
    }
    return;
}

性能测试

# Elapsed time is 0.062636 seconds.
tic; integralChannels(ones(1000, 1000, 9)); toc;
# Elapsed time is 0.085479 seconds.
tic; vl_imintegral(ones(1000, 1000, 9)); toc;
# Elapsed time is 0.128849 seconds.
tic; cumsum2D(ones(1000, 1000, 9)); toc;