基于VTK平面切割功能

medimage

// 切割介绍
// 对于一个模型的切割需要怎么办呢，想想切西瓜就知道，首先需要有一个模型、然后有一个切割平面
// 接着对于每个切割操作来更新模型，这样就可以得到切割的效果了
#include "vtkPlanes.h"
#include "vtkProperty.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkVolume.h"
#include "vtkVolumeProperty.h"
#include "vtkXMLImageDataReader.h"

#include "vtkContourFilter.h"
#include "vtkSmartPointer.h"
#include "vtkPolyDataNormals.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkOutlineFilter.h"
#include "vtkStripper.h"

#include "vtkSmoothPolyDataFilter.h"

#include <vtkSphereSource.h>
#include <vtkImagePlaneWidget.h>
#include <vtkInteractorStyleTrackballActor.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include "vtkActor.h"
#include "vtkImageFlip.h"
#include "vtkImageResample.h"
#include "vtkImageViewer.h"
#include "vtkConeSource.h"
#include "vtkBoxWidget.h"
#include "vtkTransform.h"

#include "VTKReBuild.h"
#include <vtkLineSource.h>
#include <vtkDataSetMapper.h>
#include <vector>
#include <OpenCV243.h>

#include "vtkPolyDataWriter.h"
#include "vtkPolyDataReader.h"

using namespace std;
using namespace cv;

class BuildVTKWidgetCall : public vtkCommand
{
public:

static BuildVTKWidgetCall *New()
{
  return new BuildVTKWidgetCall;
}
public:
virtual void Execute(vtkObject *caller, unsigned long eventId, void *callData)
{
     vtkImplicitPlaneWidget *pWidget = reinterpret_cast<vtkImplicitPlaneWidget*>(caller);
     if (pWidget)
     {
         vtkSmartPointer<vtkPlane> planeNew = vtkPlane::New();
         pWidget->GetPlane(planeNew);

         cliper->SetClipFunction(planeNew);
         cliper->Update();

         vtkSmartPointer<vtkPolyData> clipedData = vtkPolyData::New();
         clipedData->DeepCopy(cliper->GetOutput());

         vtkSmartPointer<vtkPolyDataMapper> coneMapper = vtkPolyDataMapper::New();
         coneMapper->SetInput(clipedData);
         coneMapper->ScalarVisibilityOff();   
         actor->SetMapper(coneMapper);
     }
}
void setCliper(vtkSmartPointer<vtkClipPolyData> other){cliper = other;}
void setPlane(vtkSmartPointer<vtkPlane> other){pPlane = other;}
void setActor(vtkSmartPointer<vtkActor> other){actor = other;}
private:
     vtkSmartPointer<vtkPlane> pPlane;
     vtkSmartPointer<vtkActor> actor;
     vtkSmartPointer<vtkClipPolyData> cliper;
};

void build3DView()
{
    vtkSmartPointer<vtkRenderer> aRenderer =
        vtkSmartPointer<vtkRenderer>::New();
    vtkSmartPointer<vtkRenderWindow> renWin =
        vtkSmartPointer<vtkRenderWindow>::New();
    renWin->AddRenderer(aRenderer);

    vtkSmartPointer<vtkRenderWindowInteractor> iren =
        vtkSmartPointer<vtkRenderWindowInteractor>::New();
    iren->SetRenderWindow(renWin);

    vtkSmartPointer<vtkJPEGReader> dicomReader =
        vtkSmartPointer<vtkJPEGReader>::New();  

    dicomReader->SetFilePrefix("C:/Users/DawnWind/Desktop/000/");
    dicomReader->SetFilePattern("%s%d.jpg");
    dicomReader->SetDataByteOrderToLittleEndian();
    dicomReader->SetDataSpacing(1, 1, 1.4);
    dicomReader->SetFileNameSliceSpacing(1);
    dicomReader->SetDataExtent(0, 209, 0, 209, 0, 29);
    dicomReader->Update();  

    vtkSmartPointer<vtkContourFilter> skinExtractor =
        vtkSmartPointer<vtkContourFilter>::New();
    skinExtractor->SetInputConnection(dicomReader->GetOutputPort());
    skinExtractor->SetValue(0, 100);    //值越大，保留的部分越少。
#ifdef TEST   
    // 这里有用到OpenCV 如果用户没有安装OpenCV那么请将与之有关的删除   
    // 这里一定要update不然下面的getpoints之类是无法取得数据的
    jpegReader->Update()   
    skinExtractor->Update();
    auto data = skinExtractor->GetOutput();
    auto points = data->GetPoints();
    auto pSize = points->GetNumberOfPoints();
    vector<Point3d> pointsGroup;
    Mat newMat = Mat::zeros(210, 210, CV_8UC1);
    int matStep = newMat.step;
    auto matData = newMat.data;
    Point2d center;
    for (int i = 0; i < pSize; i++)
    {
        double point[3];
        points->GetPoint(i, point);
        Point3d p1;
        p1.x = (point[0]);
        p1.y = (point[1]);
        p1.z = (point[2]);
        *(matData + (int)point[0] + (int)point[1] * matStep) = 255;
        pointsGroup.push_back(p1);
        center.x += (int)point[0];
        center.y += (int)point[1];
    }
    center.x /= pSize;
    center.y /= pSize;
    imshow("mat", newMat);
    //Mat dst0;
    //flip(newMat, dst0, 0);
    //imshow("dst0", dst0);
    //Mat dst1;
    //flip(newMat, dst1, 1);
    //imshow("dst1", dst1);
    //Mat dstn1;
    //flip(newMat, dstn1, -1);
    //imshow("dstn1", dstn1);
    //imwrite("a.jpg", newMat);
    // 图像本身是与原始图像成某轴对称因此不能在原图中找中心点
    waitKey();
#endif
    /**做平滑处理**/
    vtkSmartPointer<vtkSmoothPolyDataFilter> smooth = vtkSmoothPolyDataFilter::New();
    smooth->SetInput( skinExtractor->GetOutput());
    smooth->SetNumberOfIterations( 100 );

    //重新计算法向量
    vtkSmartPointer<vtkPolyDataNormals> skinNormals =
        vtkSmartPointer<vtkPolyDataNormals>::New();

    skinNormals->SetInputConnection(smooth->GetOutputPort());
    skinNormals->SetFeatureAngle(60.0);        

    vtkSmartPointer<vtkStripper> skinStripper =        //create triangle strips and/or poly-lines 为了更快的显示速度
        vtkSmartPointer<vtkStripper>::New();
    skinStripper->SetInputConnection(skinNormals->GetOutputPort());

    vtkSmartPointer<vtkPolyDataMapper> skinMapper =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    skinMapper->SetInputConnection(skinStripper->GetOutputPort());
    skinMapper->ScalarVisibilityOff();    //这样不会带颜色

    vtkSmartPointer<vtkActor> skin =
        vtkSmartPointer<vtkActor>::New();
    skin->SetMapper(skinMapper);

    vtkSmartPointer<vtkOutlineFilter> outlineData =
        vtkSmartPointer<vtkOutlineFilter>::New();
    outlineData->SetInputConnection(dicomReader->GetOutputPort());

    vtkSmartPointer<vtkPolyDataMapper> mapOutline =
        vtkSmartPointer<vtkPolyDataMapper>::New();
    mapOutline->SetInputConnection(outlineData->GetOutputPort());

    vtkSmartPointer<vtkActor> outline =
        vtkSmartPointer<vtkActor>::New();
    outline->SetMapper(mapOutline);
    outline->GetProperty()->SetColor(0, 0, 0);

    vtkSmartPointer<vtkCamera> aCamera =
        vtkSmartPointer<vtkCamera>::New();
    aCamera->SetViewUp (0, 0, -1);
    aCamera->SetPosition (0, 1, 0);
    aCamera->SetFocalPoint (0, 0, 0);
    aCamera->ComputeViewPlaneNormal();
    aCamera->Azimuth(30.0);
    aCamera->Elevation(30.0);
    aCamera->Dolly(1.5);  
    aRenderer->AddActor(outline);
    aRenderer->AddActor(skin);
    aRenderer->SetActiveCamera(aCamera);
    aRenderer->ResetCamera ();
    aRenderer->SetBackground(.2, .3, .4);
    aRenderer->ResetCameraClippingRange ();

    renWin->SetSize(640, 480);

    vtkSmartPointer<vtkInteractorStyleTrackballCamera> style =
    vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
    iren->SetInteractorStyle( style );
      /设置截面
      vtkSmartPointer<vtkClipPolyData> cliper = vtkClipPolyData::New();
      cliper->SetInput(skinStripper->GetOutput());

      // 此平面box可以通过右键来进行放大缩小处理（只有当鼠标控制区域只有切割体才单一有效）
      vtkSmartPointer<vtkImplicitPlaneWidget> implicitPlaneWidget = vtkImplicitPlaneWidget::New();
      implicitPlaneWidget->SetInteractor(iren);
      implicitPlaneWidget->SetPlaceFactor(1.25);
      //initially position the widget
      implicitPlaneWidget->SetInput(skinStripper->GetOutput());
      implicitPlaneWidget->PlaceWidget();

      //Render2
      vtkSmartPointer<vtkActor> coneSkinActor = vtkActor::New();
      coneSkinActor->SetMapper( skinMapper );

      vtkSmartPointer<vtkRenderer> rRenderer =
      vtkSmartPointer<vtkRenderer>::New();
      rRenderer->SetBackground( 0.2, 0.3, 0.5 );
      rRenderer->SetViewport(0.5, 0.0, 1.0, 1.0);
      rRenderer->AddActor(coneSkinActor);

      vtkSmartPointer<BuildVTKWidgetCall> pCall = BuildVTKWidgetCall::New();
      pCall->setActor(coneSkinActor);
      pCall->setCliper(cliper);

      renWin->AddRenderer(rRenderer);
      ///

    implicitPlaneWidget->AddObserver(vtkCommand::EndInteractionEvent, pCall);
    implicitPlaneWidget->On();
    // Render
    renWin->Render();
    // Initialize the event loop and then start it.
    iren->Initialize();

    iren->Start();
}