目前正在学习切片问题,关于vtkImageReslice如何设置旋转矩阵,我参考资料写了一个demo,
下面是代码
static double[] CrossProduct(double[] a, double[] b)
{
double[] c = new double[3];
c[0] = a[1] * b[2] - a[2] * b[1];
c[1] = a[2] * b[0] - a[0] * b[2];
c[2] = a[0] * b[1] - a[1] * b[0];
return c;
}
static double DotProduct(double[] a, double[] b)
{
double result;
result = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
return result;
}
static double Normalize(double[] v)
{
double result;
result = Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
return result;
}
static double[,] RotationMatrix(double angle, double[] u)
{
double norm = Normalize(u);
double[,] rotatinMatrix = new double[3, 3];
u[0] = u[0] / norm;
u[1] = u[1] / norm;
u[2] = u[2] / norm;
rotatinMatrix[0, 0] = Math.Cos(angle) + u[0] * u[0] * (1 - Math.Cos(angle));
rotatinMatrix[0, 1] = u[0] * u[1] * (1 - Math.Cos(angle)) - u[2] * Math.Sin(angle);
rotatinMatrix[0, 2] = -u[1] * Math.Sin(angle) + u[0] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[1, 0] = -u[2] * Math.Sin(angle) + u[0] * u[1] * (1 - Math.Cos(angle));
rotatinMatrix[1, 1] = Math.Cos(angle) + u[1] * u[1] * (1 - Math.Cos(angle));
rotatinMatrix[1, 2] = u[0] * Math.Sin(angle) + u[1] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 0] = u[1] * Math.Sin(angle) + u[0] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 1] = -u[0] * Math.Sin(angle) + u[1] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 2] = Math.Cos(angle) + u[2] * u[2] * (1 - Math.Cos(angle));
return rotatinMatrix;
}
static void Calculation(double[] vectorBefore, double[] vectorAfter)
{
double[] rotationAxis;
double rotationAngle;
double[,] rotationMatrix;
rotationAxis = CrossProduct(vectorBefore, vectorAfter);
rotationAngle = Math.Acos(DotProduct(vectorBefore, vectorAfter) / Normalize(vectorBefore) / Normalize(vectorAfter));
rotationMatrix = RotationMatrix(rotationAngle, rotationAxis);
}
使用该算法对比vtk官方给的例子
Examples/ImageProcessing/Cxx/ImageSlicing.cxx · master · VTK / VTK · GitLab (kitware.com)
我发现冠状面计算和例子中一致,但是矢状面不一样。
我的理解是求冠状面的话是求向量 (0,1,0)到向量(0,0,1)的旋转矩阵
求横切面的话是求向量(0,0,1)到向量(0,0,1)的旋转矩阵
求矢状面的话是求向量(1,0,0)到向量(0,0,1)的旋转矩阵
不知道是理解错了还是算法错了,求大佬指点。
下面是代码
static double[] CrossProduct(double[] a, double[] b)
{
double[] c = new double[3];
c[0] = a[1] * b[2] - a[2] * b[1];
c[1] = a[2] * b[0] - a[0] * b[2];
c[2] = a[0] * b[1] - a[1] * b[0];
return c;
}
static double DotProduct(double[] a, double[] b)
{
double result;
result = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
return result;
}
static double Normalize(double[] v)
{
double result;
result = Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
return result;
}
static double[,] RotationMatrix(double angle, double[] u)
{
double norm = Normalize(u);
double[,] rotatinMatrix = new double[3, 3];
u[0] = u[0] / norm;
u[1] = u[1] / norm;
u[2] = u[2] / norm;
rotatinMatrix[0, 0] = Math.Cos(angle) + u[0] * u[0] * (1 - Math.Cos(angle));
rotatinMatrix[0, 1] = u[0] * u[1] * (1 - Math.Cos(angle)) - u[2] * Math.Sin(angle);
rotatinMatrix[0, 2] = -u[1] * Math.Sin(angle) + u[0] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[1, 0] = -u[2] * Math.Sin(angle) + u[0] * u[1] * (1 - Math.Cos(angle));
rotatinMatrix[1, 1] = Math.Cos(angle) + u[1] * u[1] * (1 - Math.Cos(angle));
rotatinMatrix[1, 2] = u[0] * Math.Sin(angle) + u[1] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 0] = u[1] * Math.Sin(angle) + u[0] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 1] = -u[0] * Math.Sin(angle) + u[1] * u[2] * (1 - Math.Cos(angle));
rotatinMatrix[2, 2] = Math.Cos(angle) + u[2] * u[2] * (1 - Math.Cos(angle));
return rotatinMatrix;
}
static void Calculation(double[] vectorBefore, double[] vectorAfter)
{
double[] rotationAxis;
double rotationAngle;
double[,] rotationMatrix;
rotationAxis = CrossProduct(vectorBefore, vectorAfter);
rotationAngle = Math.Acos(DotProduct(vectorBefore, vectorAfter) / Normalize(vectorBefore) / Normalize(vectorAfter));
rotationMatrix = RotationMatrix(rotationAngle, rotationAxis);
}
使用该算法对比vtk官方给的例子
Examples/ImageProcessing/Cxx/ImageSlicing.cxx · master · VTK / VTK · GitLab (kitware.com)
我发现冠状面计算和例子中一致,但是矢状面不一样。
我的理解是求冠状面的话是求向量 (0,1,0)到向量(0,0,1)的旋转矩阵
求横切面的话是求向量(0,0,1)到向量(0,0,1)的旋转矩阵
求矢状面的话是求向量(1,0,0)到向量(0,0,1)的旋转矩阵
不知道是理解错了还是算法错了,求大佬指点。