Very long time to deform FE mesh with Thin-plate Spline Transformation

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.

  1. Extract the boundary surface points of the liver region from the whole volumetric mesh (box + liver) and take it as source points.

  2. Apply deformation to the liver surface/boundary by simply applying translation, scaling, etc, and then extract the deformed surface points as target points.

  3. Next, initialize thin-plate spline transform with this source and target point sets.

  4. 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();

This is a common problem for transforming FE meshes based on landmarks. You can use ScatteredTransform module in 3D Slicer to create a bspline transform from TPS transform, which makes transformations magnitudes faster. It comes with MIT license, so if you want then you can extract the logic of the module and create a VTK filter from it.

this actually takes 1 second to run:

from vedo import *
tetm = TetMesh('volumetric_liver_box.vtk')
geom = tetm.threshold(below=1).tomesh(fill=False).c('red').alpha(1)
source = geom.clone().clean(0.1)
target = source.clone().scale([1.4,0.7,1.4]).addPos(0.2,0.1,0.3)
arrows = Arrows(source, target, s=0.3)
deformed_geom = geom.clone().thinPlateSpline(source, target).alpha(0.2)
print(type(deformed_geom.transform))
show(geom, arrows, deformed_geom, axes=7)

<class 'vtkCommonTransformsPython.vtkThinPlateSplineTransform'>

I think the trick is to use vtkCleanPolyData with a tolerance to greatly reduce the number of landmark points.

1 Like