-
功能描述（straightened CPR）：
如图1所示，输入的是原始冠脉CTA影像，图2提供的是原始冠脉影像与提供的中心线的空间位置关系，现在需要对原始冠脉CTA影像进行拉直曲面重构，并实现重构后的实时的切面显示交互操作。
-
解决方法：
1）首先对输入的中心线进行重采样，通过调用vtk的vtkSplineFilter类进行重采样操作，
然后用vmtk的vtkvmtkCenterlineAttributesFilter和vtkvmtkCenterlineGeometry类对中心线添加后续重构需要的属性“ParallelTransportNormals”和“FrenetTangent”；
2）根据中心线上的点生成各个点在给定的冠脉影像上的切面，可以通过vtk中的vtkImageReslice类，设定好切面的大小等参数，沿着重采样后的新的中心线在每个点上生成一张切面图；
3）将各个点的切面正确地堆叠在一起生成一个完整的重构模型，可以通过自己定义一个vtkImageData对象，大小可以从前面的切面一些参数取得，然后在每次生成一个切面的时候就遍历这个切面的所有像素值，并把值赋给我们定义的新的vtkImageData对象（注：这里也尝试过布尔运算等一些其他方法，均无法实现最后的效果）；
4）添加一种交互操作能够实时地显示整条血管各个角度的切面，这里我们用vtk首先生成了一个平面，设定好平面的朝向和位置，然后通过修改鼠标滚轮事件的调用函数来实现平面的旋转，因为我们最终重构好的模型是沿着Z轴堆叠的，所以这里我们的平面也只要沿着Z轴旋转就可以了。接下去的步骤就只要把平面旋转后相对于一开始变化的角度传出来，就可以计算出平面旋转后的法向量，这样就可以更新平面的朝向，当然与模型相交的部分也同时实现了更新。
-
具体实现代码：
void Straight() {
vtkSmartPointer<vtkPolyData> line;
VtkUtil::ReadVTP(line, PATH + "/L9.vtp");
int resolution = 600;
double step_length = line->GetLength() / resolution;
vtkSmartPointer<vtkSplineFilter> spline =
vtkSmartPointer<vtkSplineFilter>::New();
spline->SetInputData(line);
spline->SetSubdivideToLength();
spline->SetLength(step_length);
spline->Update();
line = spline->GetOutput();
vtkSmartPointer<vtkvmtkCenterlineAttributesFilter> attributesFilter =
vtkSmartPointer<vtkvmtkCenterlineAttributesFilter>::New();
attributesFilter->SetInputData(line);
attributesFilter->SetAbscissasArrayName("Abscissas");
attributesFilter->SetParallelTransportNormalsArrayName("ParallelTransportNormals");
attributesFilter->Update();
vtkSmartPointer<vtkvmtkCenterlineGeometry> geometry =
vtkSmartPointer<vtkvmtkCenterlineGeometry>::New();
geometry->SetInputData(attributesFilter->GetOutput());
geometry->SetLengthArrayName("Length");
geometry->SetCurvatureArrayName("Curvature");
geometry->SetTorsionArrayName("Torsion");
geometry->SetTortuosityArrayName("Tortuosity");
geometry->SetFrenetTangentArrayName("FrenetTangent");
geometry->SetFrenetNormalArrayName("FrenetNormal");
geometry->SetFrenetBinormalArrayName("FrenetBinormal");
geometry->SetLineSmoothing(0);
geometry->SetOutputSmoothedLines(0);
geometry->SetNumberOfSmoothingIterations(100);
geometry->SetSmoothingFactor(0.1);
geometry->Update();
vtkSmartPointer<vtkImageData> Image = vtkSmartPointer<vtkImageData>::New();
Image = image_data_;
vtkSmartPointer<vtkPolyData> centerline;
centerline = geometry->GetOutput();
vtkSmartPointer<vtkImageData> out_image;
double a[3] = {0.1, 0.1, 0.1};
int size_one = 50;
_stepAngle = 45.0 / (size_one / 2);
int b[2] = {size_one, size_one};
this->CurvedMPRImages(Image, centerline, 0.0, a, b, step_length, out_image);
vtkSmartPointer<vtkXMLImageDataReader> imageSource =
vtkSmartPointer<vtkXMLImageDataReader>::New();
QString cur_path = PATH + "/final.vti";
imageSource->SetFileName(cur_path.toLocal8Bit().data());
imageSource->Update();
Image = imageSource->GetOutput();
this->_image = Image;
double bounds[6];
Image->GetBounds(bounds);
_planeWidth = bounds[5];
_planeHeight = bounds[3];
vtkSmartPointer<vtkDataSetMapper> mapper =
vtkSmartPointer<vtkDataSetMapper>::New();
mapper->SetInputData(Image);
vtkSmartPointer<vtkActor> actor =
vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
this->_finalActor = actor;
image_seeder_->WidgetsOff();
this->vmtk_renderer_->GetRenderer()->AddActor(actor);
this->vmtk_renderer_->Render();
this->vmtk_renderer_->ResetCamera();
}
void CurvedMPRImages(vtkSmartPointer<vtkImageData> &image,
vtkSmartPointer<vtkPolyData> ¢erline,
const double reslicingBackgroundLevel,
const double outputSpacing[3],
int outputSize[2],
double stepLength,
vtkSmartPointer<vtkImageData> &outImage) {
vtkDataArray *frenetTangentArray = centerline->GetPointData()->GetArray("FrenetTangent");
vtkDataArray *parallelTransportNormalsArray = centerline->GetPointData()
->GetArray("ParallelTransportNormals");
vtkPoints *linePoints = centerline->GetPoints();
int numberOfLinePoints = linePoints->GetNumberOfPoints();
int OutputExtent[6];
OutputExtent[0] = 0;
OutputExtent[1] = outputSize[0] - 1;
OutputExtent[2] = 0;
OutputExtent[3] = outputSize[1] - 1;
OutputExtent[4] = 0;
OutputExtent[5] = numberOfLinePoints - 1;
double OutputOrigin[3];
double OutputSpacing[3];
double firstPoint[3], secondPoint[3], zspace;
linePoints->GetPoint(0, firstPoint);
linePoints->GetPoint(1, secondPoint);
zspace = sqrt(vtkMath::Distance2BetweenPoints(firstPoint, secondPoint));
OutputSpacing[0] = outputSpacing[0];
OutputSpacing[1] = outputSpacing[1];
OutputSpacing[2] = zspace;
OutputOrigin[0] = -1 * ((OutputExtent[1] - OutputExtent[0]) / 2) * OutputSpacing[0];
OutputOrigin[1] = -1 * ((OutputExtent[3] - OutputExtent[2]) / 2) * OutputSpacing[1];
OutputOrigin[2] = 0;
auto reslice = vtkSmartPointer<vtkImageReslice>::New();
reslice->SetOutputDimensionality(2);
reslice->SetInputData(image);
reslice->SetInterpolationModeToCubic();
// turn off transformation of the input spacin, origin and extent, so we can define what we want
reslice->TransformInputSamplingOff();
// set the value of the voxels that are out of the input data
reslice->SetBackgroundLevel(reslicingBackgroundLevel);
// set the outputspacing
//reslice->SetOutputSpacing(this->OutputSpacing[0],this->OutputSpacing[1],this->OutputSpacing[2]);
reslice->SetOutputSpacing(outputSpacing);
// going to MPR
vtkSmartPointer<vtkImageData> final_image
= vtkSmartPointer<vtkImageData>::New();
final_image->SetDimensions(outputSize[0], outputSize[1], numberOfLinePoints);
final_image->AllocateScalars(VTK_DOUBLE, 1);
final_image->SetSpacing(outputSpacing[0], outputSpacing[1], stepLength);
int *dims = final_image->GetDimensions();
// std::cout << "Dims: " << " x: " << dims[0] << " y: " << dims[1] << " z: " << dims[2] << std::endl;
// std::cout << "Number of points: " << final_image->GetNumberOfPoints() << std::endl;
// std::cout << "Number of cells: " << final_image->GetNumberOfCells() << std::endl;
for (int slice = 0; slice < numberOfLinePoints; slice++) {
// for each slice (or point on the line)
double centerSlice[3];
linePoints->GetPoint(slice, centerSlice);
// To calculate the outputorigin & the necessarry transform
// the vectors are retreived from the array's
// t is the vector in the direction of the Centerline, so along the z-axis in the MPR volume
double t[3];
frenetTangentArray->GetTuple(slice, t);
// p is a normal of the Centerline, directed to the 'North' direction of the inputvolume,in the MPR
// volume this will be along the x-axis
double p[3];
parallelTransportNormalsArray->GetTuple(slice, p);
double tp[3];
// tp is the crossproduct of t and p, and will be directed to the 'West' direction of the inputvolume,
// in the MPR volume this will be along the y-axis
tp[0] = (t[1] * p[2] - t[2] * p[1]);
tp[1] = (t[2] * p[0] - t[0] * p[2]);
tp[2] = (t[0] * p[1] - t[1] * p[0]);
// set the axis of the first slice according to the vectors of the first point in the line
reslice->SetResliceAxesDirectionCosines(p[0], p[1], p[2], tp[0], tp[1], tp[2], t[0], t[1], t[2]);
// the firstPoint on the Centerline coresponds to the origin of the output axes unit vectors
reslice->SetResliceAxesOrigin(centerSlice[0], centerSlice[1], centerSlice[2]);
// the outputextent will be one slice
reslice->SetOutputExtent(OutputExtent[0], OutputExtent[1], OutputExtent[2], OutputExtent[3], 0, 0);
// the outputorigin in the input coordinate system is set to zero
reslice->SetOutputOrigin(OutputOrigin);
reslice->Update();
for (int y = 0; y < dims[1]; y++) {
for (int x = 0; x < dims[0]; x++) {
double *pixel = static_cast<double *>(final_image->GetScalarPointer(x, y, slice));
double value = reslice->GetOutput()->GetScalarComponentAsDouble(x, y, 0, 0);
//std::cout << "value: " << value << std::endl;
pixel[0] = value;
}
}
}
QString path = PATH + "/final.vti";
vtkSmartPointer<vtkXMLImageDataWriter> writer =
vtkSmartPointer<vtkXMLImageDataWriter>::New();
writer->SetInputData(final_image);
writer->SetFileName(path.toStdString().c_str());
writer->Write();
outImage = final_image;
}
- 补充：
1）这里读取的可以是DICOM影像序列，也可以是一个完整的vti文件；
2）此处的平面交互代码较为基础就不记录了，只要根据旋转的角度更新平面对应的法向量即可。