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;