I have created a conformal mesh in gmsh by wrapping a volumetric liver mesh with a box with two physical entities, one is to hold the entire volume and the other is to hold the liver which is inside of the box. When I save the mesh in VTK format the physical tags (1 and 2) are saved in the CELL_DATA field. This whole mesh consists of 22882 points and 135111 number of cells (i.e. tetrahedrons)
Now I have applied Thin-plate Spline transformation to the whole mesh as below steps.
-
Extract the boundary surface points of the liver region from the whole volumetric mesh (box + liver) and take it as source points.
-
Apply deformation to the liver surface/boundary by simply applying translation, scaling, etc, and then extract the deformed surface points as target points.
-
Next, initialize thin-plate spline transform with this source and target point sets.
-
Apply this TPS transformation to the whole volumetric mesh.
But it will take a very long time to produce the result (still running during the past 4 hours). Any reason for this or solution to this?
Please find the code below for your further reference. I have herewith attached the .vtk file as well.
volumetric_liver_box.zip (1.6 MB)
std::string ref_volumetric_mesh_path = "./volumetric_liver_box.vtk";
vtkSmartPointer<vtkUnstructuredGridReader> reader =
vtkSmartPointer<vtkUnstructuredGridReader>::New();
reader->SetFileName(ref_volumetric_mesh_path.c_str());
reader->ReadAllVectorsOn();
reader->ReadAllScalarsOn();
reader->Update();
// Extract ellipsoid mesh
vtkSmartPointer<vtkThreshold> threshold_ellipsoid =
vtkSmartPointer<vtkThreshold>::New();
threshold_ellipsoid->SetInputData(reader->GetOutput());
threshold_ellipsoid->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, vtkDataSetAttributes::SCALARS);
threshold_ellipsoid->ThresholdByLower(1);
threshold_ellipsoid->Update();
//geometry filter to extract the boundary
vtkSmartPointer<vtkGeometryFilter> ellipsoidGeometryFilter =
vtkSmartPointer<vtkGeometryFilter>::New();
ellipsoidGeometryFilter->SetInputData(threshold_ellipsoid->GetOutput());
ellipsoidGeometryFilter->Update();
vtkSmartPointer<vtkPolyData> polyDataForEllipsoid = ellipsoidGeometryFilter->GetOutput();
vtkSmartPointer<vtkPoints> sourcePoints = polyDataForEllipsoid->GetPoints();
// apply deformation to the surface points
vtkSmartPointer<vtkTransform> transform1 =
vtkSmartPointer<vtkTransform>::New();
transform1->Translate(0.2, 0.1, 0.3);
transform1->Scale(1.4, 0.7, 1.4);
vtkSmartPointer<vtkTransformFilter> transformFilter =
vtkSmartPointer<vtkTransformFilter>::New();
transformFilter->SetInputConnection(ellipsoidGeometryFilter->GetOutputPort());
transformFilter->SetTransform(transform1);
transformFilter->Update();
//extract deformed surface points
vtkSmartPointer<vtkPoints> targetPoints = transformFilter->GetOutput()->GetPoints();
vtkSmartPointer<vtkThinPlateSplineTransform> splineTransform =
vtkSmartPointer< vtkThinPlateSplineTransform >::New();
splineTransform->SetSourceLandmarks(sourcePoints);
splineTransform->SetTargetLandmarks(targetPoints);
splineTransform->SetBasisToR2LogR();
//apply trasnformation to the whole grid
vtkSmartPointer<vtkTransformFilter> wholeMeshTransformFilter =
vtkSmartPointer<vtkTransformFilter>::New();
wholeMeshTransformFilter->SetInputConnection(reader->GetOutputPort());
wholeMeshTransformFilter->SetTransform(splineTransform);
std::string deformed_volumetric_mesh_path = "./deformed_volumetric_liver_box.vtk";
vtkSmartPointer<vtkUnstructuredGridWriter> unstructuredGridWriter =
vtkSmartPointer<vtkUnstructuredGridWriter>::New();
unstructuredGridWriter->SetFileName(deformed_volumetric_mesh_path.c_str());
unstructuredGridWriter->SetInputConnection(wholeMeshTransformFilter->GetOutputPort());
unstructuredGridWriter->Write();