opencv 三维矩阵的转置

背景

opencv矩阵类为Mat，本文主要讨论3维矩阵，更高维度的操作方式是类似的。
三维Mat按照其header的实现分为两种：二维多通道以及三维单通道。这里之所以说是按照header来分，是因为从数据存储来看，其实没有差别（重要的是数据的维度次序），就是一组内存数据而已（当然数据类型要一样，这也是来自header信息），关键是header的信息决定了能够使用哪些操作方法。这两种不同类别，即使你的内存数据次序是一样的，由于header不同，有些操作也不能使用，比如copyTo方法。而只能使用更为低级的内存拷贝操作。

二维多通道矩阵

这种类型的Mat主要是opencv用来实现图像处理的。因此其数据组织以及header信息是针对二维图像的，比如对RGB3通道的图像。opencv提供的一些方法（和图像处理相关）也是能用于这类Mat，典型的如t()转置方法就是用来实现图像的宽高转置；再比如flip()方法就是用来实现图像宽高平面的翻转(而不涉及通道)，因此也只能用于这类Mat。
再举个最简单的例子。Mat默认的std::cout输出就只能用于这种Mat，而对于三维Mat，输出会报错。

int dims[] = {2,3,4};
cv::Mat mat(3, dims, CV_32F, cv::Scalar(0)); //全0的三维矩阵
std::cout << mat << std::endl;

运行时错误：

Error: Assertion failed (m.dims <= 2) in FormattedImpl

二位多通道Mat有rows和cols属性，分别表示图像的高和宽，其dims属性值为2。而对于三维单通道，rows和cols取值为-1。

三维单通道矩阵

既然二维多通道矩阵是opencv用来实现图像处理的主要数据结构，并且很多方法的实现也是针对这类矩阵的，而从本质上来说，它其实也是三维的（高、宽、通道）。那为什么还需要三维单通道矩阵呢？或者说，我要处理三维数据，是选择二维多通道矩阵还是选择三维单通道矩阵呢？
其实这个问题是取决于你的三维数据是什么以及你将如何使用这各个维度的数据。
首先，你要处理的就是普通的图像数据（HWC或WHC)，那么肯定二维多通道矩阵更为方便，毕竟有那么多现成的方法可以使用；
其次，即使你要处理的三维数据并非图像，但是其中的两个维度是并列关系（如平面的点，甚至是各个样本），主要的reshape变换或者转置这类变换仅限于这两个维度；而最后一个维度是需要进行运算的数据，除了对它切片和连接，它不会和其他维度进行次序的交换。那么也可以使用二维多通道这种Mat（甚至二维单通道Mat）。
最后一种情况，三维数据在不同场景下需要转置为合适的次序，从而能够便于实现所需要的高效运算，这是用二维多通道Mat就不合适了。典型的场景就是高光谱数据处理。当需要针对各个通道平面进行平面空间运算的时候，希望的数据维度是CHW或CWH；当需要需要每个平面点的各个通道数据变换运算时，HWC或WHC则更为便利。
因此，需要实现三维单通道的转置运算，而opencv并未提供这种操作方法

三维单通道矩阵转置

我们希望实现类似于numpy的transpose操作。
原3D矩阵的维度依次为(0，1，2)；则通过指定新的维度次序，比如(1，2，0)，来实现相应的转置。

实现方式一：reshape+t

对于原维度次序(0，1，2)， 如果新的维度次序为(1，2，0)或(2，0，1)；也就是说，假设我们把0，1，2看成是顺时针排列，转置结果仍然是顺时针排列，这样的话，可以使用Mat的reshape方法加转置方法t来实现3D矩阵的转置

reshape操作和二维转置测试

reshape方法并未改变数据内存，即reshape的结果赋值给一个新的Mat并不会导致数据的拷贝，而仅仅改变了header信息，该Mat和原Mat的数据内存共享

int dims[] = {2,3,4};
// 创建3维矩阵，形状为2*3*4
cv::Mat x = (cv::Mat_<uint8_t>(3,dims)<< 1, 2, 3, 4, 5, 6, 
                                        7, 8, 9, 10, 11, 12, 
                                        13, 14, 15, 16, 17, 18,
                                        19, 20, 21, 22, 23, 24);
// 对3维矩阵进行reshape，结果形状为6*4
cv:Mat y = x.reshape(1, 6);
print_mat(y, "y");
// 进行2D转置，结果形状为4*6
cv::Mat z = y.t();
print_mat(z, "z");

输出结果：

y (6 x 4)
[ 1, 2, 3, 4;
5, 6, 7, 8;
9, 10, 11, 12;
13, 14, 15, 16;
17, 18, 19, 20;
21, 22, 23, 24]
z (4 x 6)
[ 1, 5, 9, 13, 17, 21;
2, 6, 10, 14, 18, 22;
3, 7, 11, 15, 19, 23;
4, 8, 12, 16, 20, 24]

查看x和y的数据内存地址，是一样的，为0x1154f00。

转置会造成数据顺序的改变，因此需要新申请内存，因此转置的结果矩阵的数据内存地址和x、y都不同。

转置前的内存数据（16进制），可以看到顺序为1，2，3，4，…， 24：

转置后的内存数据，可以看到和打印结果一致（1，5，9，13…）：

3D转置实现

以（0，1，2）转置为（1，2，0）为例。
实现步骤：

1. 3D变2D：将（0，1，2）reshape为（0，1*2）；
2. 进行2D转置：将（0，12）变为（12，0）；
3. 2D转3D：将（1*2，0）reshape为（1，2，0），得到结果。

// 3个维度：c h w
int c=2;
int h = 3;
int w = 4;
// 创建3维矩阵，形状为2*3*4
int dims[] = {c,h,w};
cv::Mat src = (cv::Mat_<uint8_t>(3,dims)<< 1, 2, 3, 4, 5, 6, 
                                        7, 8, 9, 10, 11, 12, 
                                        13, 14, 15, 16, 17, 18,
                                        19, 20, 21, 22, 23, 24);
int c_hw[] = { c, h * w };
// 3D->2D; 2D 转置
cv::Mat dst = src.reshape(1, 2, c_hw).t();
// 2D -> 3D
int dst_dims[] = {h ,w, c };
dst= dst.reshape(1, 3, dst_dims);

实现方式二：通用方式，内存数据坐标转换

虽然上述代码可以实现3D Mat转置，但是无法从(0，1，2)转置为（1，0，2），究其原因，是因为reshape操作不能改变内存数据，因此需要使用其他更为通用的实现方式。
想象一个3D的立方体，对于其上任意一点，所谓的转置，其实就是其坐标点在在不同的坐标轴次序下如何表示的问题。因此，我们可以将每个点从原坐标次序的内存位置映射到新坐标次序要求的内存位置即可。
所以，首先我们要了解3D Mat数据在内存中是如何布局的。

Mat数据内存布局

以3D为例，Mat的数据首地址为data属性的值，矩阵任意元素M(i,j,k)的地址addr为：
$$addr(M[i,j,k])=M.data+step(0)\ast i+step(1)\ast j+k+step(2)\ast k$$
这里：
$step(0)$的意思是第1维的维度之间间隔的字节数。它等于第2维元素数量$\times$第3维的元素数量$\times$每个元素的字节数。以上述代码为例，$step(0)=3\ast 4\ast 1$
$step(1)$的意思是第2维的维度之间间隔的字节数。它等于第3维元素数量$\times$每个元素的字节数。以上述代码为例，$step(1)=4\ast 1$
$step(2)$的意思是第3维的维度之间间隔的字节数。它直接等于每个元素的字节数。以上述代码为例，$step(2)=1$

因此，我们可以这样实现转置：

1. 申明新矩阵，其大小和原矩阵一样；
2. 根据新维度次序，遍历新矩阵的每个元素（假设坐标表示为i，j，k）
3. 根据维度变换规则，计算（i，j，k）在原矩阵维度次序下面的坐标(i’, j’, k’)
4. 根据内存地址公式，根据(i’, j’, k’)得到原矩阵中的元素值，赋值到新矩阵（i，j，k）位置。

代码实现示例

• 并行取代遍历：利用opencv Mat的forEach方法，对每一个元素实现相同的操作
• at方法赋值：可以直接使用上述M.data计算内存地址进行赋值；或者利用at方法进行赋值

out_mat = Mat::zeros(3, out_ndims, CV_8U);
out_mat.forEach<int8_t>(
       [&in_mat, old_order](int8_t& val, const int* position) {
         int8_t old_v = *reinterpret_cast<int8_t*>(in_mat.data + in_mat.step[0] * position[old_order[0]] + in_mat.step[1] * position[old_order[1]] + in_mat.step[2] * position[old_order[2]]);
         val = old_v;
             });

代码以CV_8U数据类型为例，in_mat为原矩阵，out_mat为转置后的矩阵（其维度次序为out_dims数组）；old_order为新矩阵维度转换为原矩阵维度的次序数组。