OpenCV Mat的实现

由于Mat在OpneCV里的地位非常重要，这篇文章打算好好研究一番...估计时间跨度会非常长...慢慢更新...

关于Mat的定义，参看我的这篇博客 http://blog.csdn.net/gauss_acm/article/details/50808753 

先上常用的函数，因为很多函数内部都可能用到这些函数，所以这样方便看下面的复杂函数

看这些前，建议先看看我上面的博客，对Mat有一个更好的理解，这样下面看起来就轻松了

//! returns true iff the matrix data is continuous// (i.e. when there are no gaps between successive rows).// similar to CV_IS_MAT_CONT(cvmat->type)bool isContinuous() const;//! returns true if the matrix is a submatrix of another matrixbool isSubmatrix() const;//! returns element size in bytes,// similar to CV_ELEM_SIZE(cvmat->type)size_t elemSize() const;//! returns the size of element channel in bytes.size_t elemSize1() const;//! returns element type, similar to CV_MAT_TYPE(cvmat->type)int type() const;//! returns element type, similar to CV_MAT_DEPTH(cvmat->type)int depth() const;//! returns element type, similar to CV_MAT_CN(cvmat->type)int channels() const;//! returns step/elemSize1()size_t step1(int i=0) const;//! returns true if matrix data is NULLbool empty() const;//! returns the total number of matrix elementssize_t total() const;

inline bool Mat::isContinuous() const { return (flags & CONTINUOUS_FLAG) != 0; }

直接看连续标志就可以了

inline bool Mat::isSubmatrix() const { return (flags & SUBMATRIX_FLAG) != 0; }

直接看子矩阵标志就可以了，这个在ROI里很常用

inline size_t Mat::elemSize() const { return dims > 0 ? step.p[dims-1] : 0; }

返回矩阵每个元素的字节数，depth对应的数据类型*通道数，比如CV_8UC3，depth是CV_8U，是1个字节，通道数是3，所以elemSize为1*3，而这个值刚好step.p[dims-1],因为最后一维寻址刚好跳跃的间距是每个元素所占的字节数

inline size_t Mat::elemSize1() const { return CV_ELEM_SIZE1(flags); }

这个就是elemSize/通道数，也就是只考虑depth对应字节数就可以了

/* Size of each channel item,   0x124489 = 1000 0100 0100 0010 0010 0001 0001 ~ array of sizeof(arr_type_elem) */#define CV_ELEM_SIZE1(type) \    ((((sizeof(size_t)<<28)|0x8442211) >> CV_MAT_DEPTH(type)*4) & 15)/* 0x3a50 = 11 10 10 01 01 00 00 ~ array of log2(sizeof(arr_type_elem)) */#define CV_ELEM_SIZE(type) \    (CV_MAT_CN(type) << ((((sizeof(size_t)/4+1)*16384|0x3a50) >> CV_MAT_DEPTH(type)*2) & 3))

对应两个变态的宏，其实也不变态了，例如第 2个宏里的0x3a50是7种类型的编码，对于uchar和schar是单字节的，我们需要搞个1出来，而CV_8U,CV_8S是0和1，我们乘2得到0和2，0x3a50右移0位和2位，然后&3，&3相当于提取末两位，这样右移0位和2位得到的末2位都是0，所以我们把depth左移0位相当于乘1...于是就用这个方法凑出每种类型对应的末2位，这样得到了0x3a50这个数了...前面16384*2实际上就是1<<15,这个就是子矩阵的标志，或上去在这个宏里不会影响末2位的，也就是这里是多余的，可以删掉2*16384...

inline int Mat::type() const { return CV_MAT_TYPE(flags); }inline int Mat::depth() const { return CV_MAT_DEPTH(flags); }inline int Mat::channels() const { return CV_MAT_CN(flags); }

这些看 http://blog.csdn.net/gauss_acm/article/details/50785649 ，分别取type(12位二进制数)，depth(3位二进制数)和channels(数字)

inline size_t Mat::step1(int i) const { return step.p[i]/elemSize1(); }

这个感觉没什么意义...

inline size_t Mat::total() const{    if( dims <= 2 )        return rows*cols;    size_t p = 1;    for( int i = 0; i < dims; i++ )        p *= size[i];    return p;}

矩阵元素个数

inline bool Mat::empty() const { return data == 0 || total() == 0; }

矩阵为空的条件是data=null或者元素个数为0

先来看一个重要的函数

//! sets every matrix element to sMat& operator = (const Scalar& s);

这个函数规定了Mat的通道数小于等于4,所以Scalar是4维的够了，这个函数的作用就是将Mat的每个元素用Scalar赋值，如果Mat是k通道,k<=4,那么Scalar的前k维分别赋值给Mat的k维通道...注意一下，Scalar实际上是Scalar_<double>,对于Mat元素不是double型的，则系统自动使用saturate_cast进行转换,关于saturate_cast，参考我的博客 http://blog.csdn.net/gauss_acm/article/details/50811783 

具体实现

Mat& Mat::operator = (const Scalar& s){    const Mat* arrays[] = { this };    uchar* ptr;    NAryMatIterator it(arrays, &ptr, 1);    size_t size = it.size*elemSize(); //获得超平面字节数    if( s[0] == 0 && s[1] == 0 && s[2] == 0 && s[3] == 0 )    {        //全0的情况直接按字节填充0        for( size_t i = 0; i < it.nplanes; i++, ++it )            memset( ptr, 0, size );    }    else    {        if( it.nplanes > 0 ) //先填充第一个超平面        {            double scalar[12];            scalarToRawData(s, scalar, type(), 12);            size_t blockSize = 12*elemSize1();            for( size_t j = 0; j < size; j += blockSize )            {                size_t sz = MIN(blockSize, size - j);                memcpy( ptr + j, scalar, sz );            }        }        for( size_t i = 1; i < it.nplanes; i++ )        {            //填充完第一个超平面后data已经有数据了,因为ptr是指针            ++it;            memcpy( ptr, data, size ); //用第一个超平面拷贝到其他超平面        }    }    return *this;}

里面有NAryMatIterator这个类，实际上这个是一个迭代器，推荐看一下我的另一篇 http://blog.csdn.net/gauss_acm/article/details/51195572 ，看完之后，你会发现在opencv里考虑了矩阵不连续存储的情况，所以这个迭代器将会以超平面的形式遍历矩阵

我们看看scalarToRawData这个函数

void scalarToRawData(const Scalar& s, void* _buf, int type, int unroll_to){    int i, depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);    CV_Assert(cn <= 4);    switch(depth)    {    case CV_8U:        {        uchar* buf = (uchar*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<uchar>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_8S:        {        schar* buf = (schar*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<schar>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_16U:        {        ushort* buf = (ushort*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<ushort>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_16S:        {        short* buf = (short*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<short>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_32S:        {        int* buf = (int*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<int>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_32F:        {        float* buf = (float*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<float>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        }        break;    case CV_64F:        {        double* buf = (double*)_buf;        for(i = 0; i < cn; i++)            buf[i] = saturate_cast<double>(s.val[i]);        for(; i < unroll_to; i++)            buf[i] = buf[i-cn];        break;        }    default:        CV_Error(CV_StsUnsupportedFormat,"");    }}

对不同的depth采用了不同的数据转换，因为unroll_to传进来的是12，也就是每次赋值12个，这里为什么是12？因为通道数只能是1,2,3,4,所以为了不把1个元素的通道分割开来，我们需要用1,2,3,4的公倍数，而12就是最小的公倍数，举个例子，如果通道数可以是5，那么假设超平面有3个元素，那么就是3*5=15个，而此时用12去复制一次，还剩3个，然后size_t sz = MIN(blockSize, size - j);这句代码会去sz=min(12,3)=3，于是再用12的前3个去复制剩余的3个，那么我们想，这样会是正确的吗？显然不对啊，因为12中的前10个刚好是2*5，对应两个矩阵的元素，12中的后两个对应第3个元素的前两个通道，所以剩余3个不能用12的前3个来填，这个是衔接不上的...所以需要公倍数

这个复制语句我们将在很多Mat的构造函数看到...

最重要的3个create函数

//! allocates new matrix data unless the matrix already has specified size and type.// previous data is unreferenced if needed.void create(int _rows, int _cols, int _type);void create(Size _size, int _type);void create(int _ndims, const int* _sizes, int _type);

具体实现

inline void Mat::create(int _rows, int _cols, int _type){    _type &= TYPE_MASK;    //如果矩阵参数相同并且数据区域!=null不创建矩阵    if( dims <= 2 && rows == _rows && cols == _cols && type() == _type && data )        return;    int sz[] = {_rows, _cols};    create(2, sz, _type);}inline void Mat::create(Size _sz, int _type){    create(_sz.height, _sz.width, _type);}void Mat::create(int d, const int* _sizes, int _type){    int i;    CV_Assert(0 <= d && _sizes && d <= CV_MAX_DIM && _sizes);    _type = CV_MAT_TYPE(_type);    if( data && (d == dims || (d == 1 && dims <= 2)) && _type == type() )    {        if( d == 2 && rows == _sizes[0] && cols == _sizes[1] )            return;        for( i = 0; i < d; i++ )            if( size[i] != _sizes[i] )                break;        if( i == d && (d > 1 || size[1] == 1))            return;    }    release();    if( d == 0 )        return;    flags = (_type & CV_MAT_TYPE_MASK) | MAGIC_VAL;    setSize(*this, d, _sizes, 0, allocator == 0);    if( total() > 0 )    {        if( !allocator )        {            size_t total = alignSize(step.p[0]*size.p[0], (int)sizeof(*refcount));            data = datastart = (uchar*)fastMalloc(total + (int)sizeof(*refcount));            refcount = (int*)(data + total);            *refcount = 1;        }        else        {            allocator->allocate(dims, size, _type, refcount, datastart, data, step.p);            CV_Assert( step[dims-1] == (size_t)CV_ELEM_SIZE(flags) );        }    }    finalizeHdr(*this);}

其实前两个都会调用第3个，所以看第3个就可以了...

当allocator==NULL,则会分配data空间，datastart=data，引用计数设置为1... 其他函数详见下面分析

CV_MAX_DIM=32，也就是说最多创建32维的矩阵，上来先判断需要创建的矩阵是否已存在，否则release

inline void Mat::release(){    if( refcount && CV_XADD(refcount, -1) == 1 )        deallocate();    data = datastart = dataend = datalimit = 0;    size.p[0] = 0;    refcount = 0;}

void Mat::deallocate(){    if( allocator )        allocator->deallocate(refcount, datastart, data);    else    {        CV_DbgAssert(refcount != 0);        fastFree(datastart);    }}

CV_XADD可以看成后置自增运算.,可以参考我的博客 http://blog.csdn.net/gauss_acm/article/details/50971503 

也就是当refcount为1，马上减为0，就需要释放掉，调用dellocate

先看一下setSize

static inline void setSize( Mat& m, int _dims, const int* _sz,                            const size_t* _steps, bool autoSteps=false ){    CV_Assert( 0 <= _dims && _dims <= CV_MAX_DIM );    if( m.dims != _dims )    {        if( m.step.p != m.step.buf )        {            fastFree(m.step.p);            m.step.p = m.step.buf;            m.size.p = &m.rows;        }        if( _dims > 2 )        {            m.step.p = (size_t*)fastMalloc(_dims*sizeof(m.step.p[0]) + (_dims+1)*sizeof(m.size.p[0]));            m.size.p = (int*)(m.step.p + _dims) + 1;            m.size.p[-1] = _dims;            m.rows = m.cols = -1;        }    }    m.dims = _dims;    if( !_sz )        return;    size_t esz = CV_ELEM_SIZE(m.flags), total = esz;    int i;    for( i = _dims-1; i >= 0; i-- )    {        int s = _sz[i];        CV_Assert( s >= 0 );        m.size.p[i] = s;        if( _steps )            m.step.p[i] = i < _dims-1 ? _steps[i] : esz;        else if( autoSteps )        {            m.step.p[i] = total;            int64 total1 = (int64)total*s;            if( (uint64)total1 != (size_t)total1 )                CV_Error( CV_StsOutOfRange, "The total matrix size does not fit to \"size_t\" type" );            total = (size_t)total1;        }    }    if( _dims == 1 )    {        m.dims = 2;        m.cols = 1;        m.step[1] = esz;    }}

我们在读代码时，发现Mat是没有1维的情况，函数自动将1维转成2维，而第二维长度为1,其次对于2维的情况，size.p和step.p是不需要申请空间的，size.p[0]和rows的地址相同，size.p[1]与clos的地址 相同，同理step.p和buf的地址相同，对于大于2维的情况rows=cols=-1，此时buf没用了 ,size.p和step.p独立申请空间，可以发现size.p多申请了4个字节的空间，size.p[-1]用来存dims(维数)，step数组可以传进来，当然也可以不传进来，通过size数组来算，这个刚好可以用于第一次创建矩阵，autoStep可以传true进来...create里使用了allocator==0显然第一次创建allocator必然为0，所以相当于传了true...

里面还有fastMalloc和fastFree两个函数，这是opencv库分配和释放内存的函数，其中使用了指针对齐的技术，

详见我的博客 http://blog.csdn.net/gauss_acm/article/details/50971503 

create里的alignSize

/*!  Aligns buffer size by the certain number of bytes  This small inline function aligns a buffer size by the certian number of bytes by enlarging it.*/static inline size_t alignSize(size_t sz, int n){    return (sz + n-1) & -n;}

注意当n为2的幂才有意义，例如n=4,二进制100，-n的二进制位111...1100,&-n相当于截取高位为n的倍数的部分，这里当且仅当n为2的幂次，实际上就是求大于等于n的最小的n的倍数...

static void finalizeHdr(Mat& m){    updateContinuityFlag(m);    int d = m.dims;    if( d > 2 )        m.rows = m.cols = -1;    if( m.data )    {        m.datalimit = m.datastart + m.size[0]*m.step[0];        if( m.size[0] > 0 )        {            m.dataend = m.data + m.size[d-1]*m.step[d-1];            for( int i = 0; i < d-1; i++ )                m.dataend += (m.size[i] - 1)*m.step[i];        }        else            m.dataend = m.datalimit;    }    else        m.dataend = m.datalimit = 0;}

这个函数做了最后的一些成员变量的更新，可以看到当d>2时，rows=cols=-1,datastart是数据的起始位置，其值等于data，datalimit是数据区域的末地址(包括矩阵不连续存储的末尾的一些空白字节)，dataend是真正的数据的结束位置，它和datalimit的区别是当矩阵不连续存储，末尾部分空白字节不算进去了...所以只有当矩阵连续存储时，这两个值是一样的...

static void updateContinuityFlag(Mat& m){    int i, j;    for( i = 0; i < m.dims; i++ )    {        if( m.size[i] > 1 )            break;    }    for( j = m.dims-1; j > i; j-- )    {        if( m.step[j]*m.size[j] < m.step[j-1] )            break;    }    int64 t = (int64)m.step[0]*m.size[0];    if( j <= i && t == (int)t )        m.flags |= Mat::CONTINUOUS_FLAG;    else        m.flags &= ~Mat::CONTINUOUS_FLAG;}

更新连续标志，非常简单，首先找到最开始的一维，长度大于1，因为长度为1，只有一个超平面，是不会有间隔之说的(也就是肯定连续的)，然后从d-1维枚举到i，如果满足不等式则说明中间有填补的空白字符，是不连续的，于是更新flags标记...否则将flags的CONTINUOUS_FLAG(1<<14)这位清空，使用位运算就是&~CONTINUOUS_FLAG

接下来是Mat一些构造函数

//! default constructorMat();//! constructs 2D matrix of the specified size and type// (_type is CV_8UC1, CV_64FC3, CV_32SC(12) etc.)Mat(int _rows, int _cols, int _type);Mat(Size _size, int _type);//! constucts 2D matrix and fills it with the specified value _s.Mat(int _rows, int _cols, int _type, const Scalar& _s);Mat(Size _size, int _type, const Scalar& _s);//! constructs n-dimensional matrixMat(int _ndims, const int* _sizes, int _type);Mat(int _ndims, const int* _sizes, int _type, const Scalar& _s);

没什么好说的，基本上都是在内部调用create函数，Scalar则使用我上面讲的Mat = Scalar函数

inline Mat::Mat()    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){}inline Mat::Mat(int _rows, int _cols, int _type)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create(_rows, _cols, _type);}inline Mat::Mat(int _rows, int _cols, int _type, const Scalar& _s)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create(_rows, _cols, _type);    *this = _s;}inline Mat::Mat(Size _sz, int _type)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create( _sz.height, _sz.width, _type );}    inline Mat::Mat(Size _sz, int _type, const Scalar& _s)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create(_sz.height, _sz.width, _type);    *this = _s;}    inline Mat::Mat(int _dims, const int* _sz, int _type)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create(_dims, _sz, _type);}inline Mat::Mat(int _dims, const int* _sz, int _type, const Scalar& _s)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    create(_dims, _sz, _type);    *this = _s;}    

拷贝构造

//! copy constructorMat(const Mat& m);

inline Mat::Mat(const Mat& m)    : flags(m.flags), dims(m.dims), rows(m.rows), cols(m.cols), data(m.data),    refcount(m.refcount), datastart(m.datastart), dataend(m.dataend),    datalimit(m.datalimit), allocator(m.allocator), size(&rows){    if( refcount )        CV_XADD(refcount, 1);    if( m.dims <= 2 )    {        step[0] = m.step[0]; step[1] = m.step[1];    }    else    {        dims = 0;        copySize(m);    }}

当refcount不为0时，才加1，为什么呢？因为当refcount为0时，说明原矩阵只有一个矩阵头，并没有数据区域，所以拷贝构造过来当然也没有数据区域了...

注意：拷贝构造函数值拷贝了一个矩阵头...数据区域是公用的...所以当改变拷贝矩阵的数据时，原矩阵也跟着改变...要小心使用...

当维数小于等于2的情况，step复制一下就可以了，因为size在冒号语法里已初始化，否则调用copySize

void Mat::copySize(const Mat& m){    setSize(*this, m.dims, 0, 0);    for( int i = 0; i < dims; i++ )    {        size[i] = m.size[i];        step[i] = m.step[i];    }}

在copySize里调用了setSize，注意调用前dims清0，这样在setSize里step和size会申请空间...最后复制一下数据就可以了

利用用户自己构造的数据来创建Mat

//! constructor for matrix headers pointing to user-allocated dataMat(int _rows, int _cols, int _type, void* _data, size_t _step=AUTO_STEP);Mat(Size _size, int _type, void* _data, size_t _step=AUTO_STEP);Mat(int _ndims, const int* _sizes, int _type, void* _data, const size_t* _steps=0);

具体实现

inline Mat::Mat(int _rows, int _cols, int _type, void* _data, size_t _step)    : flags(MAGIC_VAL + (_type & TYPE_MASK)), dims(2), rows(_rows), cols(_cols),    data((uchar*)_data), refcount(0), datastart((uchar*)_data), dataend(0),    datalimit(0), allocator(0), size(&rows){    size_t esz = CV_ELEM_SIZE(_type), minstep = cols*esz;    if( _step == AUTO_STEP )    {        _step = minstep;        flags |= CONTINUOUS_FLAG;    }    else    {        if( rows == 1 ) _step = minstep;        CV_DbgAssert( _step >= minstep );        flags |= _step == minstep ? CONTINUOUS_FLAG : 0;    }    step[0] = _step; step[1] = esz;    datalimit = datastart + _step*rows;    dataend = datalimit - _step + minstep;}

这个函数创建的是二维矩阵,1~4通道，函数参数里data是用户创建的数据区域首地址，_step是每一行元素所占的字节数，_step默认是0，此时函数会自动根据cols*elemSize来计算...但也可以传入一个具体的值，这个值可以大于cols*elemSize，也就是用户自己分配的数据区域同样可以不连续，那么更新flag的连续性标志就可以根据这个大于来做...dataend是数据结束的地址，不包括填充的空白字节，所以减去_step，最后一行单独算，于是加上cols *elemSize，由于这块 数据区域是用户自己分配的，所以refcount不需要设置为1，也就是这块区域无需Mat自己释放...

inline Mat::Mat(Size _sz, int _type, void* _data, size_t _step)    : flags(MAGIC_VAL + (_type & TYPE_MASK)), dims(2), rows(_sz.height), cols(_sz.width),    data((uchar*)_data), refcount(0), datastart((uchar*)_data), dataend(0),    datalimit(0), allocator(0), size(&rows){    size_t esz = CV_ELEM_SIZE(_type), minstep = cols*esz;    if( _step == AUTO_STEP )    {        _step = minstep;        flags |= CONTINUOUS_FLAG;    }    else    {        if( rows == 1 ) _step = minstep;        CV_DbgAssert( _step >= minstep );        flags |= _step == minstep ? CONTINUOUS_FLAG : 0;    }    step[0] = _step; step[1] = esz;    datalimit = datastart + _step*rows;    dataend = datalimit - _step + minstep;}

第二个类似，不讲了...

Mat::Mat(int _dims, const int* _sizes, int _type, void* _data, const size_t* _steps)    : flags(MAGIC_VAL|CV_MAT_TYPE(_type)), dims(0),    rows(0), cols(0), data((uchar*)_data), refcount(0),    datastart((uchar*)_data), dataend((uchar*)_data), datalimit((uchar*)_data),    allocator(0), size(&rows){    setSize(*this, _dims, _sizes, _steps, true);    finalizeHdr(*this);}

调用setSize和finalizeHdr，注意最后一维必须要连续，因为setSize会把step[dims-1]设为elemSize的，官方文档也有说明...

inline Mat& Mat::operator = (const Mat& m){    if( this != &m )    {        if( m.refcount )            CV_XADD(m.refcount, 1);        release();        flags = m.flags;        if( dims <= 2 && m.dims <= 2 )        {            dims = m.dims;            rows = m.rows;            cols = m.cols;            step[0] = m.step[0];            step[1] = m.step[1];        }        else            copySize(m);        data = m.data;        datastart = m.datastart;        dataend = m.dataend;        datalimit = m.datalimit;        refcount = m.refcount;        allocator = m.allocator;    }    return *this;}

复制构造函数，如果两个对象指针相同，直接返回，否则先维护引用数，再维护其他域，如果小于等2维，直接拷贝，大于2为需要copySize，这个函数上面已分析...

接下来我们考虑如何截取子矩阵，同时更新flags|=SUBMATRIX_FLAG,我们不需要重新复制数据区域，只需要增加引用计数，并且我们需要修改data的首地址为子矩阵的第一个元素的地址，同时需要更新每一维的size，但是我们不需要更新step，为什么呢？因为我们引用的是之前的矩阵的数据区域，这个数据区域可能比子矩阵大，所以子矩阵不一定连续存储，所以我们计算子矩阵某个元素的地址，还是需要按照之前的矩阵的step来计算，同时我们考虑更新CONTINUOUS_FLAG标志...

//! creates a matrix header for a part of the bigger matrixMat(const Mat& m, const Range& rowRange, const Range& colRange=Range::all());Mat(const Mat& m, const Rect& roi);Mat(const Mat& m, const Range* ranges);

还有三个括号运算符

//! extracts a rectangular sub-matrix// (this is a generalized form of row, rowRange etc.)Mat operator()( Range rowRange, Range colRange ) const;Mat operator()( const Rect& roi ) const;Mat operator()( const Range* ranges ) const;

因为这之间调用关系复杂，所以我还是按照调用关系来吧

Mat::Mat(const Mat& m, const Range* ranges)    : flags(m.flags), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    int i, d = m.dims;    CV_Assert(ranges);    for( i = 0; i < d; i++ )    {        //检查Range的范围是否正确        Range r = ranges[i];        CV_Assert( r == Range::all() || (0 <= r.start && r.start < r.end && r.end <= m.size[i]) );    }    *this = m; //在复制构造函数里,引用数自动增加1    for( i = 0; i < d; i++ )    {        Range r = ranges[i];        if( r != Range::all() && r != Range(0, size.p[i]))        {            size.p[i] = r.end - r.start; //更新size的每一维            data += r.start*step.p[i]; //计算data为子矩阵的数据起始地址            flags |= SUBMATRIX_FLAG; //更新SUBMATRIX_FLAG标志        }    }    //更新连续性标志    updateContinuityFlag(*this);}

传入Range*构造n维子矩阵调用这个函数

inline Mat Mat::operator()(const Range* ranges) const{    return Mat(*this, ranges); //注意这里会有一个临时对象返回，但之后会被析构，因此引用数不会变多}

这个调用上面一个

Mat::Mat(const Mat& m, const Range& rowRange, const Range& colRange)    : flags(0), dims(0), rows(0), cols(0), data(0), refcount(0),    datastart(0), dataend(0), datalimit(0), allocator(0), size(&rows){    CV_Assert( m.dims >= 2 );    if( m.dims > 2 )    {        //大于2维,从第3维开始默认全部是Range::all()        AutoBuffer<Range> rs(m.dims);        rs[0] = rowRange;        rs[1] = colRange;        for( int i = 2; i < m.dims; i++ )            rs[i] = Range::all();        *this = m(rs); //从处rs会隐式调用转换函数变成Range*        //然后就是调用m的(Range*)运算符        return;    }    //二维情况    *this = m;    if( rowRange != Range::all() && rowRange != Range(0,rows) )    {        CV_Assert( 0 <= rowRange.start && rowRange.start <= rowRange.end && rowRange.end <= m.rows );        rows = rowRange.size();        data += step*rowRange.start; //step隐式调用转换函数,返回buf[0]        flags |= SUBMATRIX_FLAG; //更新子矩阵标记    }    if( colRange != Range::all() && colRange != Range(0,cols) )    {        CV_Assert( 0 <= colRange.start && colRange.start <= colRange.end && colRange.end <= m.cols );        cols = colRange.size();        data += colRange.start*elemSize();        //如果列数小于原矩阵,则不连续        flags &= cols < m.cols ? ~CONTINUOUS_FLAG : -1;        flags |= SUBMATRIX_FLAG; //更新子矩阵标记    }    if( rows == 1 ) //一行必然是连续的        flags |= CONTINUOUS_FLAG;    if( rows <= 0 || cols <= 0 )    {        release();        rows = cols = 0;    }}

这个函数也可用于多维，会调用上面的函数，关于AutoBuffer实际上是一个动态的数组,具体参考我的这篇博客 http://blog.csdn.net/gauss_acm/article/details/50969539

inline Mat Mat::operator()( Range rowRange, Range colRange ) const{    return Mat(*this, rowRange, colRange);}

 //相当于Range(roi.x,roi.x+roi.width),Range(roi.y,roi.y+roi.height) Mat::Mat(const Mat& m, const Rect& roi)    : flags(m.flags), dims(2), rows(roi.height), cols(roi.width),    data(m.data + roi.y*m.step[0]), refcount(m.refcount),    datastart(m.datastart), dataend(m.dataend), datalimit(m.datalimit),    allocator(m.allocator), size(&rows){    CV_Assert( m.dims <= 2 );    //子矩阵的cols小于原矩阵,不连续    flags &= roi.width < m.cols ? ~CONTINUOUS_FLAG : -1;    //一行必定连续    flags |= roi.height == 1 ? CONTINUOUS_FLAG : 0;    size_t esz = CV_ELEM_SIZE(flags);    data += roi.x*esz;    CV_Assert( 0 <= roi.x && 0 <= roi.width && roi.x + roi.width <= m.cols &&              0 <= roi.y && 0 <= roi.height && roi.y + roi.height <= m.rows );    if( refcount )        CV_XADD(refcount, 1);    if( roi.width < m.cols || roi.height < m.rows )        flags |= SUBMATRIX_FLAG;    step[0] = m.step[0]; step[1] = esz;    if( rows <= 0 || cols <= 0 )    {        release();        rows = cols = 0;    }}

inline Mat Mat::operator()( const Rect& roi ) const{ return Mat(*this, roi); }

这里做一下说明，关于子矩阵标志，当且仅当子矩阵和原矩阵不完全相同，即至少存在一维长度不同，否则不算子矩阵。

然后关于子矩阵寻址做一下严格的说明，以二维矩阵为例，设原矩阵的step信息：step[0]和step[1]，子矩阵为[x0,y0]->[x1,y1],那么对于x0<=x<=x1,y0<=y<=y1,我们可以使用y*step[0]+x*step[1]定位，但这样每一维需要记录两个端点，于是我们来做一个平移变换，以子矩阵的左上角为参考系，那么左上角地址为D=y0*step[0]+x0*step[1],于是对于[x,y],我们可以用D+(y-y0)*step[0]+(x-x0)*step[1],于是我们只需要将[x0,x1]区间映射到[0,x1-x0],[y0,y1]映射到[0,y1-y0],这样[x,y]对应了[x-x0,y-y0],这样我们每一维我们只需要记录右端点x1-x0,y1-y0,再把子矩阵左上角地址赋值给data，这样就可以正确寻址了...大概是这么一个思想...

持续更新中...