利用vtk实现管状模型沿中心线切割平面

需求就是一个管状的模型（两个开口），实现提取其中心线并沿着中心线切线方向切割模型提取平面（中心线每个点负法线方向是新平面y轴，中心线每个点切向量是新平面z轴）。用vtkvmtkPolyDataCenterlines实现中心线提取，用vtkFrenetSerretFrame实现中心校切向量和负法向量计算，用vtkCutter实现平面裁切。

1. 原始模型

2. 封闭并压缩模型（压缩为了快速提取中心线）

3. 模型中心线

4. 根据中心线法向量开始裁切

5. 所有裁切后平面点画在一起

6. 抽取其中一层比较，红色是从模型上裁切出来的。

7 实现代码

void ModelCutting();int main(int, char *[]) {ModelCutting();return 0;
}//----------------------------------------------------------
/*** @brief GetDotDistance 求两点距离* @param p* @param x* @return*/
double GetDotDistance(const double p[3], const double x[3]) {return   sqrt((p[0] - x[0]) * (p[0] - x[0]) +(p[1] - x[1]) * (p[1] - x[1]) +(p[2] - x[2]) * (p[2] - x[2]));
}//----------------------------------------------------------
/*** @brief AffineTransformation 空间坐标系变换* @param coordinate_x* @param coordinate_y* @param coordinate_z* @param center* @param point* @param outpoint*/
void AffineTransformation(const QList<double> coordinate_x,const QList<double> coordinate_y,const QList<double> coordinate_z,const QList<double> center,const QList<double> point,QList<double> &outpoint) {outpoint << (point.at(0) - center.at(0)) * coordinate_x.at(0)+ (point.at(1) - center.at(1)) * coordinate_x.at(1)+ (point.at(2) - center.at(2)) * coordinate_x.at(2) ;outpoint << (point.at(0) - center.at(0)) * coordinate_y.at(0)+ (point.at(1) - center.at(1)) * coordinate_y.at(1)+ (point.at(2) - center.at(2)) * coordinate_y.at(2);outpoint << (point.at(0) - center.at(0)) * coordinate_z.at(0)+ (point.at(1) - center.at(1)) * coordinate_z.at(1)+ (point.at(2) - center.at(2)) * coordinate_z.at(2);
}//----------------------------------------------------------
/*** @brief GetVerticalVector 已知x、z求y轴向量* @param coordinate_z z轴向量* @param coordinate_x x轴向量* @param coordinate_y y轴向量*/
void GetVerticalVector(const QList<double> coordinate_z,const QList<double> coordinate_x,QList<double> &coordinate_y) {coordinate_y[0] = coordinate_x[1] * coordinate_z[2] - coordinate_x[2] * coordinate_z[1];coordinate_y[1] = coordinate_x[2] * coordinate_z[0] - coordinate_x[0] * coordinate_z[2];coordinate_y[2] = coordinate_x[0] * coordinate_z[1] - coordinate_x[1] * coordinate_z[0];
}//----------------------------------------------------------
/*** @brief GetTwoPointUnitVector 获取两点组成的单位向量* @param point_1* @param point_2* @param unit_vctor*/
void GetTwoPointUnitVector(const QList<double> &point_1,const QList<double> &point_2,QList<double> &unit_vctor) {QList<double> vctor;vctor << point_2.at(0) - point_1.at(0)<< point_2.at(1) - point_1.at(1)<< point_2.at(2) - point_1.at(2);double distance = vctor.at(0) * vctor.at(0) +vctor.at(1) * vctor.at(1) +vctor.at(2) * vctor.at(2);distance = sqrt(distance);unit_vctor.clear();unit_vctor << vctor.at(0) / distance;unit_vctor << vctor.at(1) / distance;unit_vctor << vctor.at(2) / distance;
}//----------------------------------------------------------
/*** @brief ScatterGeneratePng* 散点图转png* @param point* @param spacing* @param size*/
void ScatterGeneratePng(const QList<QList<double>> &point,const char *out_name,const char *in_name,const double spacing = 0.4296875,const qint32 zoom = 10) {cv::Mat src = cv::imread(in_name);// flip(src, src, -1); //rotate 180qint32 sisze = src.rows * zoom;double new_spacing = spacing / zoom;double center = (sisze * new_spacing) / 2.0;cv::resize(src, src, cv::Size(sisze, sisze));foreach (auto var, point) {qint32 x = sisze - static_cast<qint32>((center + var.at(0)) / new_spacing );qint32 y = static_cast<qint32>((center - var.at(1)) / new_spacing );src.at<cv::Vec3b>(y, x)[2] = 255;}cv::imwrite(out_name, src);
}//----------------------------------------------------------
/*** @brief ComputingCenterLines* 提取中心线* @param surface* @return*/
vtkSmartPointer<vtkPolyData> ComputingCenterLines(vtkSmartPointer<vtkPolyData> surface) {qDebug() << "ComputingCenterLines begin";// 模型压缩vtkNew<vtkDecimatePro> decimation1;decimation1->SetInputData(surface);decimation1->SetTargetReduction(0.95);decimation1->Update();// 表面整理vtkNew<vtkCleanPolyData> surface_cleaner;surface_cleaner->SetInputData(decimation1->GetOutput());surface_cleaner->Update();// 三角形相交检查vtkNew<vtkTriangleFilter> surface_triangulator;surface_triangulator->SetInputConnection(surface_cleaner->GetOutputPort());surface_triangulator->PassLinesOff();surface_triangulator->PassVertsOff();surface_triangulator->Update();vtkSmartPointer<vtkPolyData> centerline_input_surface =surface_triangulator->GetOutput();vtkSmartPointer<vtkIdList> cap_center_ids = nullptr;vtkSmartPointer<vtkvmtkCapPolyData> surface_capper;// 表面封闭surface_capper = vtkSmartPointer<vtkvmtkCapPolyData>::New();surface_capper->SetInputConnection(surface_triangulator->GetOutputPort());surface_capper->SetDisplacement(0);surface_capper->SetInPlaneDisplacement(0);surface_capper->Update();centerline_input_surface = surface_capper->GetOutput();cap_center_ids = surface_capper->GetCapCenterIds();VtkUtil::ShowVtkDebugPolydata(centerline_input_surface,"Closed model  " +QString::number(cap_center_ids->GetNumberOfIds()) +"  holes");// 计算中心线vtkNew<vtkIdList> inlet_seed_ids, outlet_seed_ids;inlet_seed_ids->InsertNextId(0);outlet_seed_ids->InsertNextId(1);vtkNew<vtkvmtkPolyDataCenterlines> centerline_filter;centerline_filter->SetInputData(centerline_input_surface);centerline_filter->SetCapCenterIds(cap_center_ids);centerline_filter->SetSourceSeedIds(inlet_seed_ids);centerline_filter->SetTargetSeedIds(outlet_seed_ids);centerline_filter->SetRadiusArrayName("MaximumInscribedSphereRadius");centerline_filter->SetCostFunction("1/R");centerline_filter->SetFlipNormals(0);centerline_filter->SetAppendEndPointsToCenterlines(0);// 从端点开始centerline_filter->SetSimplifyVoronoi(0);centerline_filter->SetCenterlineResampling(0);centerline_filter->SetResamplingStepLength(1.0);centerline_filter->Update();qDebug() << "ComputingCenterLines end";return centerline_filter->GetOutput();
}//----------------------------------------------------------
/*** @brief SplinePoints 平滑点* @param adjust_list* @param spline_list* @param resolution* @return*/
bool SplinePoints(const QList<QList<double>> &adjust_list,QList<QList<double>> &spline_list, qint32 resolution) {vtkNew<vtkPoints> points;for (qint32 i = 0; i < adjust_list.size(); ++i) {points->InsertPoint(static_cast<vtkIdType>(i),adjust_list[i][0], adjust_list[i][1], adjust_list[i][2]);}vtkNew<vtkCellArray> lines;lines->InsertNextCell(adjust_list.size());vtkNew<vtkPolyData> poly_data;for (qint32 i = 0; i < adjust_list.size(); ++i) {lines->InsertCellPoint(static_cast<vtkIdType>(i));}poly_data->SetPoints(points);poly_data->SetLines(lines);vtkNew<vtkCardinalSpline> spline;spline->SetLeftConstraint(2);spline->SetLeftValue(0.0);spline->SetRightConstraint(2);spline->SetRightValue(0.0);vtkNew<vtkSplineFilter> filter;filter->SetInputData(poly_data);filter->SetNumberOfSubdivisions(resolution - 1);filter->SetSpline(spline);filter->Update();double temp_point[3];for (qint32 i = 0; i < filter->GetOutput()->GetNumberOfPoints(); ++i) {filter->GetOutput()->GetPoint(i, temp_point);spline_list.append({temp_point[0], temp_point[1], temp_point[2]});}return true;
}//----------------------------------------------------------
/*** @brief ModelCutting* 模型切割*/
void ModelCutting() {QList<QList<double>> ivus_edge_points_;qint32 guide_number = 155;// 读取模型vtkSmartPointer<vtkPolyData> surface;VtkUtil::ReadSTL(surface, "/home/arteryflow/图片/DicomData/ProCTR/mesh.stl");VtkUtil::ShowVtkDebugPolydata(surface);// 提取中心线vtkSmartPointer<vtkPolyData> center_line =vtkSmartPointer<vtkPolyData>::New();center_line = ComputingCenterLines(surface);VtkUtil::ShowVtkDebugPolydata(center_line);vtkSmartPointer<vtkPoints> points_center = vtkSmartPointer<vtkPoints>::New();for(qint32 id = 0; id < center_line->GetNumberOfPoints(); id++) {points_center->InsertNextPoint(center_line->GetPoint(center_line->GetNumberOfPoints() - 1 - id)[0],center_line->GetPoint(center_line->GetNumberOfPoints() - 1 - id)[1],center_line->GetPoint(center_line->GetNumberOfPoints() - 1 - id)[2]);}vtkNew<vtkParametricSpline> spline_center;spline_center->SetPoints(points_center);vtkNew<vtkParametricFunctionSource> functionSource_center ;functionSource_center->SetParametricFunction(spline_center);functionSource_center->SetUResolution(guide_number - 1);functionSource_center->SetVResolution(guide_number - 1);functionSource_center->SetWResolution(guide_number - 1);functionSource_center->Update();center_line = functionSource_center->GetOutput();VtkUtil::ShowVtkDebugPolydata(center_line,surface,"check");//vtkNew<vtkFrenetSerretFrame> frame;frame->SetInputData(center_line);frame->ConsistentNormalsOn();frame->Update();// 副法向量 x轴方向frame->GetOutput()->GetPointData()->SetActiveVectors("FSBinormals");vtkSmartPointer<vtkDataArray> guide_auxiliary =frame->GetOutput()->GetPointData()->GetVectors();// 导丝切向量 z轴方向frame->GetOutput()->GetPointData()->SetActiveVectors("FSTangents");vtkSmartPointer<vtkDataArray> guide_tange_vector =frame->GetOutput()->GetPointData()->GetVectors();//创建切割平面vtkNew<vtkPlane> plane;vtkNew<vtkCutter> cutter;for(qint32 i = 0; i < guide_number ; ++i) {double tmpdouble[3];// 中心点center_line->GetPoint(i, tmpdouble);QList<double> center = {tmpdouble[0], tmpdouble[1], tmpdouble[2]};// 新坐标z轴 单位向量guide_tange_vector->GetTuple(i, tmpdouble);QList<double> coordinate_z = {tmpdouble[0], tmpdouble[1], tmpdouble[2]};// 新坐标x轴 单位向量guide_auxiliary->GetTuple(i, tmpdouble);QList<double> coordinate_x = {tmpdouble[0], tmpdouble[1], tmpdouble[2]};// 新坐标y轴 单位向量QList<double> coordinate_y = {0, 0, 0};GetVerticalVector(coordinate_z, coordinate_x, coordinate_y);//double p[3];center_line->GetPoint(i, p);frame->GetOutput()->GetPointData()->SetActiveVectors("FSTangents");vtkDataArray *ptNormals2 = frame->GetOutput()->GetPointData()->GetVectors();double value[3];ptNormals2->GetTuple(i, value);plane->SetOrigin(p[0], p[1], p[2]); //设置切割平面起点plane->SetNormal(value[0], value[1], value[2]); //设置切割方向为X方向cutter->SetCutFunction(plane);//设置切割平面cutter->SetInputData(surface);//设置模型cutter->GenerateCutScalarsOn();cutter->Update();vtkSmartPointer<vtkPolyData> ResultPoly =  cutter->GetOutput();qint64 n = ResultPoly->GetNumberOfPoints();QList<QList<double>> in_points;for(qint32 j = 0; j < n; j++) {double x[3];ResultPoly->GetPoint(j, x); //获取顶点坐标if(GetDotDistance(p, x) < 5) {QList<double> src_point = {x[0], x[1], x[2]};QList<double> out_points;AffineTransformation(coordinate_x, coordinate_y, coordinate_z, center, src_point, out_points);in_points << out_points;}}SplinePoints(in_points, in_points, 200);ScatterGeneratePng(in_points, QString("/home/arteryflow/桌面/tmp/%1.png").arg(guide_number - i).toLocal8Bit().data(), QString("/home/arteryflow/图片/DicomData/ProCTR//save-mutiVOI(21x21)/%1.bmp").arg(guide_number - i).toLocal8Bit().data());ivus_edge_points_ << in_points;}vtkNew<vtkPoints> points;vtkNew<vtkCellArray> vertices;for (qint32 i = 0; i < ivus_edge_points_.size(); ++i) {points->InsertNextPoint(ivus_edge_points_[i][0],ivus_edge_points_[i][1],ivus_edge_points_[i][2]);vertices->InsertNextCell(1);vertices->InsertCellPoint(i);}vtkNew<vtkPolyData> poly_data_tmp;poly_data_tmp->SetPoints(points);poly_data_tmp->SetVerts(vertices);VtkUtil::ShowVtkDebugPolydata( poly_data_tmp, "check");
}

利用vtk实现管状模型沿中心线切割平面相关推荐

利用vtk+cgal+openmesh（或者第三方格式转换软件）做牙齿模型
如何建立牙齿模型? 使用vtk+cgal+openmesh实现模型建立与修复优化 vtk实现dicom图片生成模型,运用阈值图像法生成各类模型,不同部位牙齿和骨头的分离阈值稍有变化,运用滑动条实现阈值 ...
Intel发布神经网络压缩库Distiller：快速利用前沿算法压缩PyTorch模型
Intel发布神经网络压缩库Distiller:快速利用前沿算法压缩PyTorch模型原文:https://blog.csdn.net/u011808673/article/details/8079 ...
DL之Attention-ED：基于TF NMT利用带有Attention的 ED模型训练、测试(中英文平行语料库)实现将英文翻译为中文的LSTM翻译模型过程全记录
DL之Attention-ED:基于TF NMT利用带有Attention的 ED模型训练(中英文平行语料库)实现将英文翻译为中文的LSTM翻译模型过程全记录目录测试输出结果模型监控训练过程全 ...
联邦学习【分布式机器学习技术】【①各客户端从服务器下载全局模型；②各客户端训练本地数据得到本地模型；③各客户端上传本地模型到中心服务器；④中心服务器接收各方数据后进行加权聚合操作，得全局模型】
随着计算机算力的提升,机器学习作为海量数据的分析处理技术,已经广泛服务于人类社会. 然而,机器学习技术的发展过程中面临两大挑战: 一是数据安全难以得到保障,隐私数据泄露问题亟待解决: 二是网络安全隔离 ...
python绘制3d图-python3利用Axes3D库画3D模型图
Python3利用Axes3D库画3D模型图,供大家参考,具体内容如下最近在学习机器学习相关的算法,用python实现.自己实现两个特征的线性回归,用Axes3D库进行建模. python代码 im ...
CV之IS：利用pixellib库基于deeplabv3_xception模型对《庆余年》片段实现语义分割/图像分割简单代码全实现
CV之IS:利用pixellib库基于deeplabv3_xception模型对<庆余年>片段实现语义分割/图像分割简单代码全实现目录利用pixellib库基于deeplabv3_xc ...
CV之IS：利用pixellib库基于mask_rcnn_coco模型对《庆余年》片段实现实例分割简单代码全实现
CV之IS:利用pixellib库基于mask_rcnn_coco模型对<庆余年>片段实现实例分割简单代码全实现目录利用pixellib库的instance_segmentation函 ...
ML之4PolyR：利用四次多项式回归4PolyR模型+两种正则化(Lasso/Ridge)在披萨数据集上拟合(train)、价格回归预测(test)
ML之4PolyR:利用四次多项式回归4PolyR模型+两种正则化(Lasso/Ridge)在披萨数据集上拟合(train).价格回归预测(test) 目录输出结果设计思路核心代码输出结果设 ...
EL之Bagging：kaggle比赛之利用泰坦尼克号数据集建立Bagging模型对每个人进行获救是否预测
EL之Bagging:kaggle比赛之利用泰坦尼克号数据集建立Bagging模型对每个人进行获救是否预测目录输出结果设计思路核心代码输出结果设计思路核心代码 bagging_clf = ...
TF之pix2pix：基于TF利用Facades数据集训练pix2pix模型、测试并进行生成过程全记录
TF之pix2pix:基于TF利用Facades数据集训练pix2pix模型.测试并进行生成过程全记录目录 TB监控 1.SCALARS 2.IMAGES 3.GRAPHS 4.DISTRIBUTI ...

利用vtk实现管状模型沿中心线切割平面

1. 原始模型

2. 封闭并压缩模型（压缩为了快速提取中心线）

3. 模型中心线

4. 根据中心线法向量开始裁切

5. 所有裁切后平面点画在一起

6. 抽取其中一层比较，红色是从模型上裁切出来的。

7 实现代码

利用vtk实现管状模型沿中心线切割平面相关推荐

最新文章

热门文章