VTK maximum density projection MIP

#include <iostream>
#include <vector>
#include <Eigen>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp>   
#include <vtkDICOMImageReader.h>
#include <vtkGPUVolumeRayCastMapper.h>
#include <vtkImageData.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkSmartPointer.h>
 
#include <vtkImageShiftScale.h>
#include "testVTK.h"
using namespace std;
using namespace cv;
 
void GetVector1(vtkPlaneSource* planeSource, double v1[3])
{
	double* p1 = planeSource->GetPoint1(a);double* o = planeSource->GetOrigin(a); v1[0] = p1[0] - o[0];
	v1[1] = p1[1] - o[1];
	v1[2] = p1[2] - o[2];
}
void GetVector2(vtkPlaneSource* planeSource, double v2[3])
{
	double* p2 = planeSource->GetPoint2(a);double* o = planeSource->GetOrigin(a); v2[0] = p2[0] - o[0];
	v2[1] = p2[1] - o[1];
	v2[2] = p2[2] - o[2];
}
int main( )
{
	vtkSmartPointer<vtkPlaneSource> mpSlicePlanes = vtkSmartPointer<vtkPlaneSource>::New(a); vtkSmartPointer<vtkImageViewer2> view = vtkSmartPointer<vtkImageViewer2>::New(a); vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New(a); vtkSmartPointer<vtkRenderWindow> renw = vtkSmartPointer<vtkRenderWindow>::New(a); renw->SetSize(800.800);
	vtkSmartPointer<vtkRenderWindowInteractor> interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New(a); vtkSmartPointer<vtkDICOMImageReader> dicom = vtkSmartPointer<vtkDICOMImageReader>::New(a); dicom->SetDirectoryName("D:\\software\\SLC");
	dicom->Update(a); view->SetRenderer(renderer);
	view->SetRenderWindow(renw);
	view->SetupInteractor(interactor);
	double spacing[3];/ / interval
	dicom->GetOutput() - >GetSpacing(spacing);
	double origin[3];/ / the origin
	dicom->GetOutput() - >GetOrigin(origin);
	int extent[6];/ / range
	dicom->GetOutput() - >GetExtent(extent);
	
	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]);
	double xmin = origin[0] + spacing[0] * (extent[0] - 0.5);
	double xmax = origin[0] + spacing[0] * (extent[1] + 0.5);
	double ymin = origin[1] + spacing[1] * (extent[2] - 0.5);
	double ymax = origin[1] + spacing[1] * (extent[3] + 0.5);
	double zmin = origin[2] + spacing[2] * (extent[4] - 0.5);
	double zmax = origin[2] + spacing[2] * (extent[5] + 0.5);
	double xcenter = center[0];
	double ycenter = center[1];
	double zcenter = center[2];
 
	// axial
	mpSlicePlanes->SetOrigin(xmin, ymin, zcenter);
	mpSlicePlanes->SetPoint1(xmax, ymin, zcenter);
	mpSlicePlanes->SetPoint2(xmin, ymax, zcenter);
	vtkSmartPointer<vtkImageSlabReslice> reslice = vtkSmartPointer<vtkImageSlabReslice>::New(a); reslice->SetInputConnection(dicom->GetOutputPort());
		reslice->SetBlendModeToMax(a); reslice->SetSlabThickness(1000);
		double planeAxis1[3];
		double planeAxis2[3];
		GetVector1(mpSlicePlanes, planeAxis1);
		GetVector2(mpSlicePlanes, planeAxis2);
		double planeSizeX = vtkMath::Normalize(planeAxis1);
		double planeSizeY = vtkMath::Normalize(planeAxis2);
		double normal[3];
		mpSlicePlanes->GetNormal(normal);
		vtkSmartPointer<vtkMatrix4x4> resliceAxes = vtkSmartPointer<vtkMatrix4x4>::New(a); resliceAxes->Identity(a);for (int i = 0; i < 3; i++)
		{
			resliceAxes->SetElement(0, i, planeAxis1[i]);
			resliceAxes->SetElement(1, i, planeAxis2[i]);
			resliceAxes->SetElement(2, i, normal[i]);
		}
 
		double planeOrigin[4];
		mpSlicePlanes->GetOrigin(planeOrigin);
		planeOrigin[3] = 1.0;
		resliceAxes->Transpose(a); resliceAxes->SetElement(0.3, planeOrigin[0]);
		resliceAxes->SetElement(1.3, planeOrigin[1]);
		resliceAxes->SetElement(2.3, planeOrigin[2]);
		reslice->SetResliceAxes(resliceAxes);
		reslice->Update(a); vtkSmartPointer<vtkImageData> imagedata = reslice->GetOutput(a); view->SetInputData(imagedata);
		cout << "mip";
		view->SetColorLevel(102);
		view->SetColorWindow(4095);
		view->GetRenderer() - >ResetCamera(a); interactor->SetRenderWindow(renw);
	renw->Render(a); interactor->Start(a); }Copy the code