vtk读取CT序列mip投影

最大密度投影（maximal intensity projection，MIP）

是一种应用广泛的CT及MR图像后处理技术。
MIP 运用透视法获得二维图像，即通过计算沿着被扫描物每条射线上所遇到的最大密度像素而产生的。当光线束通过一段组织的原始图像时，图像中密度最大的像素被保留，并被投影到一个二维平面上，从而形成MIP重建图像。

MIP 能反应相应像素的X线衰减值，较小的密度变化也能在MIP图像上显示，能很好地显示血管的狭窄、扩张、充盈缺损及区分血管壁上的钙化与血管腔内的对比剂。(来自百度百科)
由于MIP 显示的是一定层厚图像中CT值最高的体素，所以变化层厚也会对图像产生影响。

可以看到，MIP 图像中的血管连续性更好。

vtkImageSlabReslice

vtk中实现mip主要关注vtkImageSlabReslice这个类，此类派生自vtkImageReslice。很像vtkImageReslice，它重新切片数据。它是多线程的。它以三维图像作为输入，并沿某个方向生成二维带厚度的MPR。

mip实现需要两个步骤：

1.设置要切的方向，通过 SetResliceAxes 或 SetResliceAxesDirectionCosines 方法设置重新切片轴方向余弦。
2.设置切的厚度。
看一下这个类的源码，继承vtkImageReslice

BlendMode 切换minip，mip，max
SlabThickness 设置厚度
SlabResolution 设置厚度的采样间距
NumBlendSamplePoints 代表切片的个数，可以根据SlabThickness 和SlabResolution计算，默认为1，就是一个切片。
这三个切换状态的函数，这几个枚举代表最大投影，均值，最小，求和。
 

可以看到构造函数设置了属性的默认值

RequestInformation函数介绍只能调用父类的函数，但是还执行了一些 初始化的内容，计算了一下NumBlendSamplePoints ，根据SlabThickness 和SlabResolution。
 

最后用父类的SetResliceAxes函数设置矩阵方向。

vtk mip的实现

int main(){
    //读取CT序列
    vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();
    reader->SetDirectoryName("D:\\SLC");
    reader->Update();

    int extent[6];//维度
    reader->GetOutput()->GetExtent(extent);
    double spacing[3];//间隔
    reader->GetOutput()->GetSpacing(spacing);
    double origin[3];//原点
    reader->GetOutput()->GetOrigin(origin);
    //计算体数据中心点
    double center[3];
    center[0] = origin[0] + spacing[0] * 0.5 * (extent[0] + extent[1]);
    center[1] = origin[1] + spacing[1] * 0.5 * (extent[2] + extent[3]);
    center[2] = origin[2] + spacing[2] * 0.5 * (extent[4] + extent[5]);

    //axial 横断面的矩阵
    double m_axialMatrix[16] = {
           1, 0, 0, 0,
           0, 1, 0, 0,
           0, 0, 1, 0,
           0, 0, 0, 1
    };
    vtkNew<vtkMatrix4x4> resliceAxes;
    resliceAxes->DeepCopy(m_axialMatrix);
    resliceAxes->SetElement(0, 3, center[0]);
    resliceAxes->SetElement(1, 3, center[1]);
    resliceAxes->SetElement(2, 3, center[2]);
    //通过extent截取对应的体数据
    vtkSmartPointer<vtkExtractVOI> extractVOI = vtkSmartPointer<vtkExtractVOI>::New();
    extractVOI->SetInputConnection(reader->GetOutputPort());
    extractVOI->SetVOI(extent);
    extractVOI->Update();
    //mip
    vtkSmartPointer<vtkImageSlabReslice> reslice = vtkSmartPointer<vtkImageSlabReslice>::New();
    reslice->SetInputConnection(extractVOI->GetOutputPort());//截取的体数据
    reslice->SetOutputDimensionality(2);//设置输出为2维图片
    reslice->SetInterpolationModeToNearestNeighbor();//设置插值方式为最邻近插值
    reslice->SetBlendModeToMax();//设置为获取最大投影，其他方式还有均值和最小
    reslice->SetSlabThickness(extent[5]);//设置厚度
    reslice->SetResliceAxes(resliceAxes);//设置矩阵
    reslice->Update();
    
    vtkSmartPointer<vtkImageData> imageData = reslice->GetOutput();
    vtkSmartPointer<vtkImageFlip>flip = vtkSmartPointer<vtkImageFlip>::New();
    flip->SetInputData(imageData);
    flip->SetFilteredAxis(1);//y轴为1，x轴为0，z轴为2；
    flip->Update();
    //这里给了一个窗宽窗位和归一化，不然没法显示
    unsigned char* data = GetImageScalarData(flip->GetOutput(), 1024, 4095);
    int dims1[3];
    imageData->GetDimensions(dims1); 
    //用opencv显示了一下，用vtkImageViewer2要更方便，可以直接设置窗宽窗位
    Mat srcMat(dims1[1], dims1[0], CV_8UC1, data);

    imshow("mip", srcMat);
    waitKey(0);
    return EXIT_SUCCESS;
    }

上面的例子是用的全部的z值厚度做的mip，这个厚度可以自己设置，其实就是extent[5]。当然也可以自己变换矩阵，得到任意方向任意厚度的mip。