看这部分也是由于突发奇想,想实现一个C++的神经网络,那么就要一个多维数组类型吧,于是又要处理数组类型转换和数组运算吧,又不想使用模板(写不好),然后就想到了OpenCV,找了源码进行学习(拷贝)...,当然这个坑已经好久没填了,数组类型基本完成了,张量类型和计算就打算参考网上的自动微分教程写成数组类型的二次封装,但是我很懒Orz
分析
类型转换
在convert.simd.hpp中,定义了如下的宏和一个具体的转换实现函数
#define DEF_CVT_FUNC(suffix, cvtfunc, _Ts, _Td, _Twvec) \
static void cvt##suffix(const uchar* src_, size_t sstep, const uchar*, size_t, \
uchar* dst_, size_t dstep, Size size, void*) \
{ \
CV_INSTRUMENT_REGION(); \
const _Ts* src = (const _Ts*)src_; \
_Td* dst = (_Td*)dst_; \
cvtfunc<_Ts, _Td, _Twvec>(src, sstep, dst, dstep, size); \
}
template<typename _Ts, typename _Td, typename _Twvec> static inline void
cvt_( const _Ts* src, size_t sstep, _Td* dst, size_t dstep, Size size )
{
sstep /= sizeof(src[0]);
dstep /= sizeof(dst[0]);
for( int i = 0; i < size.height; i++, src += sstep, dst += dstep )
{
int j = 0;
#if CV_SIMD
const int VECSZ = _Twvec::nlanes*2;
for( ; j < size.width; j += VECSZ )
{
if( j > size.width - VECSZ )
{
if( j == 0 || src == (_Ts*)dst )
break;
j = size.width - VECSZ;
}
_Twvec v0, v1;
vx_load_pair_as(src + j, v0, v1);
v_store_pair_as(dst + j, v0, v1);
}
#endif
for( ; j < size.width; j++ )
dst[j] = saturate_cast<_Td>(src[j]);
}
}
这样我们就能很方便的通过宏来定义各种类型之间的转换:
////////////////////// 8u -> ... ////////////////////////
DEF_CVT_FUNC(8u8s, cvt_, uchar, schar, v_int16)
DEF_CVT_FUNC(8u16u, cvt_, uchar, ushort, v_uint16)
DEF_CVT_FUNC(8u16s, cvt_, uchar, short, v_int16)
DEF_CVT_FUNC(8u32s, cvt_, uchar, int, v_int32)
DEF_CVT_FUNC(8u32f, cvt_, uchar, float, v_float32)
DEF_CVT_FUNC(8u64f, cvt_, uchar, double, v_int32)
DEF_CVT_FUNC(8u16f, cvt1_, uchar, float16_t, v_float32)
//... Other type
最后定义一个转换函数,根据类别获取具体的转换函数(比如 CV_8U
就对应0)
BinaryFunc getConvertFunc(int sdepth, int ddepth)
{
static BinaryFunc cvtTab[][8]; //..
}
其中的 BinaryFunc
是一个函数指针
// precomp.hpp
typedef void (*BinaryFunc)(const uchar* src1, size_t step1,
const uchar* src2, size_t step2,
uchar* dst, size_t step, Size sz,
void*);
具体到如何进行转换,如果不使用SIMD,转换方法定义在了saturate.hpp,基本就是对 saturate_cast<T>
这个模板进行各种特化,以实现类型转换时的截断。
而使用SIMD的话,这方面OpenCV在实现上一层套了一层,翻了半天才看到具体的实现代码,基本就是使用 _mm256_cvtpd_ps
等运算,这种运算自带截断处理,不需要在人工添加。可能是由于不同平台等原因,才需要进行封装吧。贴一个我写的吧,其中的 alpha
和 s
是对应了OpenCV的 convertTo
函数
// 8u32f
template<>
void cvt_(const uint8 *src, float *dst, double alpha, double s, size_t size) {
#if 0
int blocks = size / 8 * 8;
int remain = size - blocks;
__m256d alpha = _mm256_set_pd(alpha_, alpha_, alpha_, alpha_);
__m256d beta = _mm256_set_pd(beta_, beta_, beta_, beta_);
for (int i = 0; i < blocks; i += 8)
{
__m128i srcData = _mm_load_si128((__m128i *)(src + i));
__m128i xlo = _mm_unpacklo_epi16(srcData, _mm_set1_epi16(0));
__m128i xhi = _mm_unpackhi_epi16(srcData, _mm_set1_epi16(0));
__m256d ylo = _mm256_cvtepi32_pd(xlo);
__m256d yhi = _mm256_cvtepi32_pd(xhi);
__m256d res1 = _mm256_add_pd(_mm256_mul_pd(ylo, alpha), beta);
__m256d res2 = _mm256_add_pd(_mm256_mul_pd(yhi, alpha), beta);
_mm_store_ps(dst + i, _mm256_cvtpd_ps(res1));
_mm_store_ps(dst + i + 4, _mm256_cvtpd_ps(res2));
}
for (int i = size - remain; i < size; ++i)
{
dst[i] = (float)(src[i] * alpha_ + beta_);
}
#else
for (size_t i = 0; i < size; ++i) {
dst[i] = saturate_cast<float>(src[i] * alpha + s);
}
#endif
}
矩阵运算
为了减少不必要的运算,OpenCV引入了 MatExpr
类型,用于记录矩阵运算,而不实际进行计算
先说说目的,例如我们要进行 A + 1 + 2
的运算,显然我们不必计算2次,而只需要计算 A + 3
就行了,这样的计算还有一些,例如 (A + 1) * 2
可以写成 2 * A + 2
, Mat::zeros(3, 3, CV_8U) + 1
简单的说,为了记录操作类型,定义了 MatOp
类,各种加减乘除都继承于它,如 MatOp_AddEx
(实际使用时是定义了static的全局变量指针,这样才能实现多态)
以加法为例,当我们进行形如 A + B
, A + 1
类型时,我们会得到一个这样的一个结构(不是完全按照OpenCV说,只是为了演示)
MaExpr e = MatExpr(Mat A, Mat B, double alphe, double beta, double s, MatOp *addOp)
表示了 alpha * A + beta * B(如果B存在的话) + s
(注:alpha
, beta
和 s
的意义依运算而改变,例如乘法就只需要 alpha
和 s
)
当需要实际运算时,就通过 e.op->assign(e, res)
来得到结果
如果一个 MatExpr
类型和 Mat
类型进行运算,将 Mat
通过 identity
运算包装一下,最终是 MatExpr
和 MatExpr
运算的形式
MatExpr
类型和 MatExpr
进行运算时,以加法为例,会进行 e1.op->Add(e1, e2, e3)
, e1
的op会来判断是否需要把结果计算出来,如果需要,就算出结果(比如 (A*B)+(C*D)
),然后得到一个新的 MatExpr
,最终得到结果则通过 e.op->assign(e, mat)
方法
(实际实现还是很麻烦的,感觉很难通过文章解释明白,自己也是看的似懂非懂的,要考虑很多情况,每个运算类型都要考虑到,特别是乘除法。有些相同的处理方法放在了基类 MatOp
中,为了支持逻辑运算等,还加入了 MatOp_Bin
,不过这个我还没看,实际我的代码里也没有加入这些运算,连乘除法都还只支持逐元素点乘XD)
具体实现可以参考 matrix_expressions.cpp
哦,还有一点,矩阵的运算和类型转换挺像的,也通过宏定义了一堆不同类型的函数,再利用函数指针进行封装,但是稍微有点不同,放一段我写的代码吧,像 arithm_add8u
就是通过宏定义的函数,内部再通过模板函数进行计算
static ArithmFunc *getAddTab() {
static ArithmFunc addTab[] = {
(arithm_add8u), (arithm_add8s), (arithm_add16u), (arithm_add16s),
(arithm_add32s), (arithm_add64s), (arithm_add32f), (arithm_add64f)
};
return addTab;
}
void check(const Array &src1, const Array &src2, Array &dst, bool bUsed2) {
if (bUsed2) {
assert(src1.equalShape(src2) && src1.type() == src2.type());
if (&dst != &src1 && &dst != &src2 && (!src1.equalShape(dst) || dst.type() != src1.type())) {
// std::cout << "dst.reallocate(src1.shape(), src1.type());" << std::endl;
dst.reallocate(src1.shape(), src1.type());
}
} else {
if (&dst != &src1 && (!src1.equalShape(dst) || dst.type() != src1.type())) {
// std::cout << "dst.reallocate(src1.shape(), src1.type());" << std::endl;
dst.reallocate(src1.shape(), src1.type());
}
}
}
void arithm_op(const Array &src1, const Array &src2, Array &dst, ArithmFunc *tab, bool bUsed2, double *var = nullptr) {
check(src1, src2, dst, bUsed2);
// dst没有data说明src1, src2都是空数组
if (!dst.data()) {
return;
}
tab[dst.type()](src1.data(), src2.data(), dst.data(), src1.size(), var);
}
// 最外层调用的是这个函数
void add(const Array &src1, const Array &src2, Array &dst) {
arithm_op(src1, src2, dst, getAddTab(), true);
}