-
功能描述（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）此處的平面交互代碼較為基礎就不記錄了，只要根據旋轉的角度更新平面對應的法向量即可。