基于VTK的MPR实现

多平面重建（MPR）是CT三维数据呈现的重要内容，其在三维数据的任一点空间位置，采用XY、XZ和YZ三个平面切空间数据得到三个切面分别为轴状面、冠状面和矢状面，而且X/Y/Z三个轴可以绕三维坐标原点任意旋转。

以下我们探讨一下利用visualization toolkit （VTK）如何来实现MRP的三视图。

vtkCubeAxesActor 单位 vtk mpr_计算机视觉


在VTK中有一个类vtkImageReslice，其堪称几何图形过滤器中的手术刀。他可以排列，旋转，翻转，缩放，重新采样，变形, 还可以在体数据中提取图像切片卷。

  1. vtkImageReslice可以通过调用其 SetResliceAxesDirectionCosines() 方法指定切片的方向，调用其SetResliceAxesOrigin()方法切片经过的原点。
  2. vtkImageReslice也可以通过调用其SetResliceAxes()方法指定切面的坐标系矩阵，格式为vtkMatrix4x4形如，
{
1, 0, 0, 10
0, 0, 1, 20
0,-1, 0, 10
0, 0, 0, 1
}

其中第一列（1，0，0）表示X方向向量坐标，第二列（0，0，-1）表示Y方向向量坐标，第三列（0，1，0）表示Z方向向量坐标，第四列（10，20，10）表示坐标系原点。另外X,Y,Z方向向量需满足满足右手法则，即X叉乘Y为Z，Y叉乘Z为X，Z叉乘X为Y。
3. 以SetResliceAxes()为例vtkImageReslice应用代码如下：

//轴状面的切面坐标系矩阵
double Axial[16] = {
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1 };
AxialResliceMatrix = vtkMatrix4x4::New();
AxialResliceMatrix->DeepCopy(Axial);
vtkSmartPointer<vtkImageReslice> myImageReslice = vtkSmartPointer<vtkImageReslice>::New();
//设置体数据来源
myImageReslice ->SetInputData(SourceData);
//设置vtkImageReslice的切面坐标系矩阵
myImageReslice ->SetResliceAxes(AxialResliceMatrix);
//设置输出是一个切片，而不是一个卷
myImageReslice ->SetOutputDimensionality(2);
//设置切面算法的插值方式为线性插值
myImageReslice ->SetInterpolationModeToLinear();

//设置颜色查找表
vtkSmartPointer<vtkLookupTable> myLookupTable = vtkSmartPointer<vtkLookupTable>::New();
myLookupTable ->SetRange(0, 1000);
myLookupTable ->SetValueRange(0.0, 1.0);
myLookupTable ->SetSaturationRange(0.0, 0.0);
myLookupTable ->SetRampToLinear();
myLookupTable ->Build();

//映射到切面图像
vtkSmartPointer<vtkImageMapToColors> myMapToColors = vtkSmartPointer<vtkImageMapToColors>::New();
myMapToColors ->SetLookupTable(myLookupTable);
myMapToColors ->SetInputConnection(myImageReslice ->GetOutputPort());
myMapToColors ->Update();

//设置连接到图像Actor，此时将myImageActor连接设置好的管道中的vtkRenderer上就完成了
vtkSmartPointer<vtkImageActor> myImageActor = vtkSmartPointer<vtkImageActor>::New();
myImageActor->SetInputData(myMapToColors ->GetOutput());
  1. 要实现三视图需要建立三个vtkImageReslice，设置三个正交切面，假如以（0，0，0）为坐标系原点则，三视图初始切面矩阵如下：
//轴状面
double Axial[16] = {
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1 };
//冠状面
double Coronal [16] = {
		1, 0, 0, 0,
		0, 0, -1, 0,
		0, 1, 0, 0,
		0, 0, 0, 1 };
//矢状面
double Sagittal [16] = {
	    0, 0, 1, 0,
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 0, 1 };
  1. 要实现三视图随三维坐标轴的移动旋转的显示，则需要调动VTK的交互机制。
    这里采用VTK中的观察者/命令模式，将所有交互都放置在一个vtkCommand交互类中。
    由于需要实现三视图的联动，这里采用将单个视图的管线收到视图类中，并把三视图的对象传递给交互类。
    视图类：
class MPRWidget : public QVTKOpenGLNativeWidget
{
	Q_OBJECT
private:
	enum Direction
	{
		X_DIRECTION = 0,
		Y_DIRECTION,
		Z_DIRECTION,
		N_DIRECTION
	};
	vtkSmartPointer<vtkMatrix4x4> AxialResliceMatrix;
	vtkSmartPointer<vtkMatrix4x4> CoronalResliceMatrix;
	vtkSmartPointer<vtkMatrix4x4> SagittalResliceMatrix;

public:
	MPRWidget(QWidget *parent);
	~MPRWidget();

private:
//当前视图的vtkWindow窗口
	vtkSmartPointer<vtkGenericOpenGLRenderWindow> myWindow;

public:
//当前视图的vtkRenderer舞台
	vtkSmartPointer<vtkRenderer> myRenderer;
//当前视图的vtkInteractor交互器
	vtkSmartPointer<vtkRenderWindowInteractor> myInteractor;
//当前视图十字光标横轴
	vtkSmartPointer<vtkActor> HorizontalLineActor;
//当前视图十字光标竖轴
	vtkSmartPointer<vtkActor> VerticalLineActor;
//当前视图切面算法类
	vtkSmartPointer<vtkImageReslice> myImageReslice;
	vtkSmartPointer<vtkImageMapToColors> myMapToColors;
	vtkSmartPointer<vtkImageActor> myImageActor;
//对应其他两个视图的对象指针
	MPRWidget* MyWidget1;
	MPRWidget* MyWidget2;

public:
//输入体数据和切面方向，其中direction为0代表矢状面，1代表冠状面，2代表轴状面
	void SetInputData(vtkImageData* data, int direction);

};

视图类：

class MPRCallBack:public vtkCommand
{
private:
	int Direction;
	MPRWidget *MyWidget;
	MPRWidget *MyWidget1;
	MPRWidget *MyWidget2;
private:
//切面轴操作判断标志，旋转
	bool translateEnable;
//切面轴操作判断标志，平移
	bool rotateEnable;
//切面轴选中判断标志，竖轴选中
	bool vActorSlected;
//切面轴选中判断标志，横轴选中
	bool hActorSlected;
public:
	MPRCallBack();
	~MPRCallBack();

	static  MPRCallBack* New()
	{
		return new MPRCallBack;
	}

	void Execute(vtkObject *caller, unsigned long event, void *callData);

//设置对应视图切面方向
	void SetDirection(int);
//传递三视图对象
	void SetMyWidget(MPRWidget*);

	void SetMyWidget1(MPRWidget*);

	void SetMyWidget2(MPRWidget*);
};
  1. 三视图联动，包括平移和旋转
    A.平移：对于平移只需要在vtkImageReslice的切面轴矩阵上右乘平移矩阵，以下为对应不同视图的平移矩阵
//冠状面视图
double translateC[16] = {
1,0,0,DisX,
0,1,0,0,
0,0,1,DisY,
0,0,0,1 };
//矢状面视图
double translateS[16] = {
1,0,0,0,
0,1,0,DisX,
0,0,1,DisY,
0,0,0,1 };
//轴状面视图
double translateA[16] = {
1,0,0,DisX,
0,1,0,DisY,
0,0,1,0,
0,0,0,1 };
vtkMatrix4x4* matrix = MyWidget-> myImageReslice ->GetResliceAxes();
vtkMatrix4x4* translateM = vtkMatrix4x4::New();
vtkMatrix4x4* resultM = vtkMatrix4x4::New();
translateM->DeepCopy(translateA);
vtkMatrix4x4::Multiply4x4(matrix, translateM, resultM);
matrix->DeepCopy(resultM);
translateM->Delete();
resultM->Delete();

B.旋转：对于旋转只需要在vtkImageReslice的切面矩阵上右乘旋转矩阵，以下为对应不同视图的旋转矩阵

//矢状面视图
double rotateX[16] = {
1,0,0,0,
0,cosD,-1 * sinD,0,
0,sinD,cosD,0,
0,0,0,1 };
//轴状面视图
double rotateZ[16] = {
cosD,-1 * sinD,0,0,
sinD,cosD,0,0,
0,0,1,0,
0,0,0,1 };
//冠状面视图
double rotateYR[16] = {
cosD,0,-1 * sinD,0,
0,1,0,0,
sinD,0,cosD,0,
0,0,0,1 };
vtkMatrix4x4* rotateM = vtkMatrix4x4::New();
vtkMatrix4x4* resultM = vtkMatrix4x4::New();
rotateM->DeepCopy(rotateX);
vtkMatrix4x4::Multiply4x4(matrix2, rotateM, resultM);
matrix2->DeepCopy(resultM);
rotateM->Delete();
resultM->Delete();
更新切面矩阵后需要采用以下方触发管线计算
MyWidget->MyMapToColors->Update();
MyWidget->MyRenderer->GetRenderWindow()->Render();

以上程序可以实现切面提取但如果要实现三视图的十字光标与切面同步，还需要计算当一个视图中光标移动后，所有切面与其他视图光标应该相应做什么样的移动或旋转。
假设轴状面光标旋转角度为a，矢状面光标旋转角度为b，冠状面光标旋转角度为c;
当轴状面光标平移X方向deltaX和Y方向deltaY，则切面坐标原点将移动
X方向 DisX = COS(a) * deltaX + SIN(a) * deltaY;
Y方向 DisY = COS(a) * deltaY - SIN(a) * deltaX;
但是矢状面中光标将平移X方向COS(a) * deltaY和Y方向SIN(a)* deltaY；冠状面中光标将平移X方向COS(a) * deltaX和Y方向SIN(a)* deltaX。
注1：当矢状面与冠状面中光标移动时对应另两个视图中光标也会变化可通过类似几何关系计算。
注2：当光标旋转时另外视图中光标有时需要发生平移，也可通过几何关系进行计算，进行相应补偿。
后续如果有时间，我将新开帖子对相关几何计算进行探讨!