diff --git a/.gitignore b/.gitignore index 68bc17f..2f19259 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,7 @@ __pycache__/ *$py.class # C extensions -*.so +build/*.so # Distribution / packaging .Python @@ -158,3 +158,4 @@ cython_debug/ # and can be added to the global gitignore or merged into this file. For a more nuclear # option (not recommended) you can uncomment the following to ignore the entire idea folder. #.idea/ +.vscode diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..dd37c90 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,27 @@ +cmake_minimum_required(VERSION 3.16.3...3.19.7 FATAL_ERROR) + +project(SlicerBoneMorphing) + +#----------------------------------------------------------------------------- +# Extension meta-information +set(EXTENSION_HOMEPAGE "https://www.slicer.org/wiki/Documentation/Nightly/Extensions/SlicerBoneMorphing") +set(EXTENSION_CATEGORY "Morphing") +set(EXTENSION_CONTRIBUTORS "Jan Heres (West Bohemian University)") +set(EXTENSION_DESCRIPTION "This extensions allows the user to generate and morph bone meshes based on its CT scans") +set(EXTENSION_ICONURL "https://www.example.com/Slicer/Extensions/SlicerBoneMorphing.png") +set(EXTENSION_SCREENSHOTURLS "https://www.example.com/Slicer/Extensions/SlicerBoneMorphing/Screenshots/1.png") +set(EXTENSION_DEPENDS "NA") # Specified as a list or "NA" if no dependencies + +#----------------------------------------------------------------------------- +# Extension dependencies +find_package(Slicer REQUIRED) +include(${Slicer_USE_FILE}) + +#----------------------------------------------------------------------------- +# Extension modules +add_subdirectory(SlicerBoneMorphing) +## NEXT_MODULE + +#----------------------------------------------------------------------------- +include(${Slicer_EXTENSION_GENERATE_CONFIG}) +include(${Slicer_EXTENSION_CPACK}) diff --git a/README.md b/README.md index b606956..cd1e5c5 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,126 @@ # SlicerBoneMorphing -Extentsion for 3D Slicer for bone mesh morphing. + +Extension for 3D Slicer for bone mesh morphing. + +At the moment, this module specializes for the *humerus* bone. + +## Special thanks <3 +Special thanks goes to my wonder colleague Eva C. Herbst (@evaherbst) and Arthur Port (@agporto) for creating the initial idea and their huge help during the development of this module! +Also, I would like to thank O. Hirose (@ohirose) for the research on BCPD/GBCPD and it's implementation (can be found [here](https://github.com/ohirose/bcpd)) + + +## Installation +**Supported platforms:** +- Linux (x86_64) +- Windows (x86_64) +- MacOS (both x86_64 and ARM; Slicer runs through Rosetta on ARM-based Macs) + + +Steps: +- Download the latest ZIP package from Releases +- Extract the ZIP contents to your desired folder +- Open up 3D Slicer, go to Edit -> Application Settings +- In the modules section, add the extracted contents' path to "Additional Module Paths" +- Restart 3D Slicer + +> **DISCLAIMER! After restarting, the installation process will begin. If there are any Python modules not available in Slicer, they will be installed, so the startup will take SIGNIFICANTLY MORE amount of time. Do not be scared, this is intended behaviour.** + +## Usage +After a successful install, the module will be available in the **Morphing** section. +When switching to the module, you should be greeted with the following UI: + +

+ +

+ +The UI consists of **4** main sections +- Input +- Preprocessing +- Generation +- Postprocessing + +## Architecture + +To be added. + +## Module sections + +### Input section +This section is self-explanatory. Here, you choose two input models: +- Source = Source for the generation; This is the model that represents the +- Target = Model which is non-complete => Needs its missing portions generated + +### Preprocessing section + +#### Point cloud preprocessing +Before the generation process, we usually want to preprocess the model. +Reasons could be to remove unwanted **outliers** or to smooth out the models. +First of all, the model is converted to a *point cloud* to be able to effectively perform some preprocessing steps. + +Then, a downsampling is performed. For this, you can configure the threshold for downsampling by the following parameter: +- **Downsampling distance threshold** + +After the downsampling, we compute the normals of the point cloud. It uses a radius for which the normals are calculated and maximum number of neighbours. This can be adjusted with the following parameter: +- **Normals estimation radius** +- **Normals estimation max neighbours** + +Also, we want calculate a *Fast point feature histogram* to somehow encode the local geometric properties. This method takes in two following parameters: +- **FPFH search radius** - Radius in which the FPFH is calculated in +- **FPFH max neighbours** - Maximum number of neighbours taken into account + +#### Registration section +At this moment, we have our point clouds preprocessed and ready for the next step, which is the "registration" section. +Here we try to define and calculate how and how much we need to adjust the source mesh to match the target one. + +For this, we will use the downsampled point clouds with their corresponding FPFHs from the previous step. +The concrete method we use is called **RANSAC**. It uses a process called "repeated random sub-sampling" to mitigate the effect of outliers and rotational differences as much as possible. +The behaviour of this algorithm can be adjusted by the following parameters: +- **Max iterations** +- **Distance threshold** - same meaning as in previous steps +- **Fitness threshold** - the lowest fitness between the point clouds to be accepted. The lower, the higher chance of finding a good fit. The higher, higher the chance that either *max iterations* are reached + +The result of the *RANSAC* algorithm is a bit "raw". To get the best possible fit, we perform the **ICP registration algorithm** upon the result. +This can be tuned by the following parameter: +- **ICP Distance threshold** - same meaning as in previous steps + +### Generation section +Since we now have a preprocessed meshes and with defined transformations from the *source* to the *target*, we can proceed to the **generation section**. + +For this purpose, we use a method called [Bayesian coherent point drift](https://github.com/ohirose/bcpd). It falls into the *non-rigid registration* category of algorithms, which actually performs the deformation of the mesh to increase the fit of the source. +It takes in both meshes and deforms the source into the target, similarly as we've already done in the [Registration section](#preprocessing-section). +Due to the problem that BCPD allows for "unrealistic" deformations, we have done the pre-registration steps, which lets us mitigate the chance of getting into unrealistic deformations. + +Now, BCPD allows for very fine adjustments of its behaviours using lots of different parameters. For the exact description of their effects, please refer to the documentation [here](https://github.com/ohirose/bcpd/blob/master/README.md). + +> **Note: You do NOT have to perform any kind of installation process, the BCPD and its geodesic variant are already pre-built and preconfigured for immediate using in this module.** + +**Not implemented options:** +- Terminal output +- File output + +### Postprocessing section +After our models have been merged successfully, we still want to apply a slight amount of postprocessing to reach the most optimal results. +We are basically using a bit of **filtering and smoothing** to the meshes. +For these, we let you modify the following parameters: +- **Clustering scaling** - Scaled size of voxel for within vertices that are clustered together (additionally refer to [here](http://www.open3d.org/docs/0.7.0/python_api/open3d.geometry.simplify_vertex_clustering.html) +- **Smoothing iterations** - Number of iterations of mesh smoothing + +After the whole process is done, both the generated mesh (source transformed into target, standalone) and the merged mesh (generated meshes merged with the target; "combined model") are import back into the current Slicer scene. + +

+ +

+ +## FAQ +To be added. + +## Troubleshooting +To be added. + +## Contributors + + + + + + diff --git a/SlicerBoneMorphing.png b/SlicerBoneMorphing.png new file mode 100644 index 0000000..6aae6ab Binary files /dev/null and b/SlicerBoneMorphing.png differ diff --git a/SlicerBoneMorphing/CMakeLists.txt b/SlicerBoneMorphing/CMakeLists.txt new file mode 100644 index 0000000..be3697d --- /dev/null +++ b/SlicerBoneMorphing/CMakeLists.txt @@ -0,0 +1,31 @@ +#----------------------------------------------------------------------------- +set(MODULE_NAME SlicerBoneMorphing) + +#----------------------------------------------------------------------------- +set(MODULE_PYTHON_SCRIPTS + ${MODULE_NAME}.py + ) + +set(MODULE_PYTHON_RESOURCES + Resources/Icons/${MODULE_NAME}.png + Resources/UI/${MODULE_NAME}.ui + ) + +#----------------------------------------------------------------------------- +slicerMacroBuildScriptedModule( + NAME ${MODULE_NAME} + SCRIPTS ${MODULE_PYTHON_SCRIPTS} + RESOURCES ${MODULE_PYTHON_RESOURCES} + WITH_GENERIC_TESTS + ) + +#----------------------------------------------------------------------------- +if(BUILD_TESTING) + + # Register the unittest subclass in the main script as a ctest. + # Note that the test will also be available at runtime. + slicer_add_python_unittest(SCRIPT ${MODULE_NAME}.py) + + # Additional build-time testing + add_subdirectory(testing) +endif() diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_linux_x86_64 b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_linux_x86_64 new file mode 100755 index 0000000..1919f0e Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_linux_x86_64 differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_arm b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_arm new file mode 100755 index 0000000..6bf0e74 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_arm differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_x86_64 b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_x86_64 new file mode 100755 index 0000000..eab931d Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_macos_x86_64 differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_win32.exe b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_win32.exe new file mode 100755 index 0000000..2ffeccd Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/bcpd_win32.exe differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/libblas.dll b/SlicerBoneMorphing/Resources/BCPD/exec/libblas.dll new file mode 100755 index 0000000..8ae667c Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/libblas.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/libgcc_s_dw2-1.dll b/SlicerBoneMorphing/Resources/BCPD/exec/libgcc_s_dw2-1.dll new file mode 100755 index 0000000..99a76cf Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/libgcc_s_dw2-1.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/libgfortran-3.dll b/SlicerBoneMorphing/Resources/BCPD/exec/libgfortran-3.dll new file mode 100755 index 0000000..5a5c455 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/libgfortran-3.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/libgomp-1.dll b/SlicerBoneMorphing/Resources/BCPD/exec/libgomp-1.dll new file mode 100755 index 0000000..56f6c5d Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/libgomp-1.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/liblapack.dll b/SlicerBoneMorphing/Resources/BCPD/exec/liblapack.dll new file mode 100755 index 0000000..65cf4af Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/liblapack.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/libquadmath-0.dll b/SlicerBoneMorphing/Resources/BCPD/exec/libquadmath-0.dll new file mode 100755 index 0000000..b812ff8 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/libquadmath-0.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/exec/pthreadGC-3.dll b/SlicerBoneMorphing/Resources/BCPD/exec/pthreadGC-3.dll new file mode 100755 index 0000000..8e1b542 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/exec/pthreadGC-3.dll differ diff --git a/SlicerBoneMorphing/Resources/BCPD/libbcpd_linux_x86_64.so b/SlicerBoneMorphing/Resources/BCPD/libbcpd_linux_x86_64.so new file mode 100755 index 0000000..23ecb23 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/libbcpd_linux_x86_64.so differ diff --git a/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_arm.dylib b/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_arm.dylib new file mode 100755 index 0000000..4dbb782 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_arm.dylib differ diff --git a/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_x86_64.dylib b/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_x86_64.dylib new file mode 100755 index 0000000..037c4b9 Binary files /dev/null and b/SlicerBoneMorphing/Resources/BCPD/libbcpd_macos_x86_64.dylib differ diff --git a/SlicerBoneMorphing/Resources/Icons/SlicerBoneMorhping.png b/SlicerBoneMorphing/Resources/Icons/SlicerBoneMorhping.png new file mode 100644 index 0000000..5d83ab4 Binary files /dev/null and b/SlicerBoneMorphing/Resources/Icons/SlicerBoneMorhping.png differ diff --git a/SlicerBoneMorphing/Resources/Samples/O35 M humerus.stl b/SlicerBoneMorphing/Resources/Samples/O35 M humerus.stl new file mode 100644 index 0000000..2a40ccc Binary files /dev/null and b/SlicerBoneMorphing/Resources/Samples/O35 M humerus.stl differ diff --git a/SlicerBoneMorphing/Resources/Samples/U35 F humerus partial.stl b/SlicerBoneMorphing/Resources/Samples/U35 F humerus partial.stl new file mode 100644 index 0000000..8aac146 Binary files /dev/null and b/SlicerBoneMorphing/Resources/Samples/U35 F humerus partial.stl differ diff --git a/SlicerBoneMorphing/Resources/UI/SlicerBoneMorphing.ui b/SlicerBoneMorphing/Resources/UI/SlicerBoneMorphing.ui new file mode 100644 index 0000000..cc309ea --- /dev/null +++ b/SlicerBoneMorphing/Resources/UI/SlicerBoneMorphing.ui @@ -0,0 +1,1341 @@ + + + UI + + + + 0 + 0 + 564 + 1257 + + + + + + + Input + + + true + + + false + + + + + + Select a node that represents the reference model + + + Source: + + + + + + + false + + + Select a node that represents the reference model + + + + vtkMRMLModelNode + + + + true + + + + + + + + + + + + + Select a node that represents the target model (the one to be deformed) + + + Target: + + + + + + + Select a node that represents the target model (the one to be deformed) + + + + vtkMRMLModelNode + + + + true + + + + + + + + + + + + + + + + Preprocessing parameters + + + false + + + true + + + + + + Point cloud preprocessing + + + + + + The higher the number, the lower the distance threshold for a neighbour point to be considered + + + Downsampling distance threshold: + + + + + + + The higher the number, the lower the distance threshold for the RANSAC registration + + + 100.000000000000000 + + + 55.000000000000000 + + + + + + + Normals estimation radius: + + + + + + + 100.000000000000000 + + + 4.000000000000000 + + + + + + + Normals estimation max neighbours: + + + + + + + 1 + + + 1000 + + + 100 + + + + + + + FPFH search radius: + + + + + + + 100.000000000000000 + + + 10.000000000000000 + + + + + + + FPFH max neighbours: + + + + + + + 1 + + + 1000 + + + 30 + + + + + + + + + + Registration + + + + + + Distance threshold: + + + + + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + The higher the number, the longer the registration process will take, but higher the chance to converge to a fit solution + + + Max iterations: + + + + + + + The higher the number, the longer the registration process will take, but higher the chance to converge to a fit solution + + + 1 + + + 30 + + + + + + + Fitness threshold: + + + + + + + 8 + + + 1.000000000000000 + + + 0.000100000000000 + + + 0.999000000000000 + + + + + + + ICP Distance threshold: + + + + + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + + + + + + + BCPD parameters + + + false + + + true + + + + + + Tuning parameters + + + true + + + false + + + + + + Omega. Outlier probability in (0,1). + + + Outlier probability (omega): + + + + + + + Omega. Outlier probability in (0,1). + + + + + + 3 + + + 0.001000000000000 + + + 0.999000000000000 + + + 0.001000000000000 + + + 0.100000000000000 + + + + + + + Lambda. Positive. It controls the expected length of deformation vectors. Smaller is longer. + + + Deformation vector length (lambda): + + + + + + + 3 + + + 100.000000000000000 + + + 10.000000000000000 + + + + + + + Beta. Positive. It controls the range where deformation vectors are smoothed. + + + Deformation smoothing range (beta): + + + + + + + 3 + + + 100.000000000000000 + + + 10.000000000000000 + + + + + + + Gamma. Positive. It defines the randomness of the point matching at the beginning of the optimization. + + + Point matching randomness (gamma): + + + + + + + Gamma. Positive. It defines the randomness of the point matching at the beginning of the optimization. + + + 3 + + + 100.000000000000000 + + + 0.100000000000000 + + + 0.100000000000000 + + + + + + + Kappa. Positive. It controls the randomness of mixing coefficients. + + + Mixing coefficient randomness (kappa): + + + + + + + Kappa. Positive. It controls the randomness of mixing coefficients. Value 1000 is considered "infinity" + + + 3 + + + 1000.000000000000000 + + + 1000.000000000000000 + + + + + + + + + + false + + + Show more granular controls over the BCPD algorithm. DISCLAIMER: USE WITH CAUTION! + + + Show advanced controls + + + + + + + Advanced controls + + + + + + Kernel parameters + + + true + + + false + + + true + + + + + + Geodesic Kernel + + + + + + Tau. The rate controlling the balance between geodesic and Gaussian kernels. + + + Kernel balance rate (tau): + + + + + + + I have an input mesh + + + + + + + The file that defines a triangle mesh. + + + Input mesh path: + + + + + + + The number of neighbors for each node, required for k-NN graph construction. + + + kNN neighbours + + + + + + + The number of neighbors for each node, required for k-NN graph construction. + + + 10 + + + + + + + The file that defines a triangle mesh. + + + Absolute path to the input mesh file + + + + + + + The radius that defines neighbors for each node, required for k-NN graph construction. + + + kNN neighbor radius + + + + + + + The radius that defines neighbors for each node, required for k-NN graph construction. + + + 3 + + + 100.000000000000000 + + + 10.000000000000000 + + + + + + + Beta. Positive. Gaussian function's width. + + + Gaussian function width (beta): + + + + + + + Beta. Positive. Gaussian function's width. + + + 3 + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + K tilde. Positive. Rank constraint on G. + + + Rank constraint (K tilde): + + + + + + + K tilde. Positive. Rank constraint on G. + + + 3 + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + Epsilon. Positive. Acceptable condition number of G. + + + Acceptable condition number (epsilon): + + + + + + + Epsilon. Positive. Acceptable condition number of G. + + + 3 + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + 100.000000000000000 + + + 1.000000000000000 + + + + + + + + + + Standard kernel + + + + + + Kernel type: + + + + + + + + + + + + + Kernel type: + + + + + + + + + + + + + Acceleration + + + true + + + false + + + + + + Mode: + + + + + + + Manual acceleration methods. For optimal convergence, use both Nystorm and KD Tree search + + + Manual acceleration + + + false + + + + + + Nystorm method + + + true + + + true + + + + + + Samples for computing G + + + + + + + 1000000 + + + 70 + + + + + + + Samples for computing J + + + + + + + 1000000 + + + 300 + + + + + + + Random number seed for the Nystrom method. Reproducibility is guaranteed if the same number is specified. + + + Randomness seed + + + + + + + Random number seed for the Nystrom method. Reproducibility is guaranteed if the same number is specified. + + + 1 + + + 1000000 + + + + + + + + + + KD Tree Search + + + true + + + true + + + + + + Scale factor of sigma that defines areas to search for neighbors. + + + Scale factor: + + + + + + + Scale factor of sigma that defines areas to search for neighbors. + + + 3 + + + 100.000000000000000 + + + 1.000000000000000 + + + 7.000000000000000 + + + + + + + Maximum radius to search for neighbors. + + + Maximum radius + + + + + + + Maximum radius to search for neighbors. + + + 0.150000000000000 + + + + + + + The value of sigma at which the KD tree search is turned on. + + + Sigma threshold: + + + + + + + The value of sigma at which the KD tree search is turned on. + + + 3 + + + 100.000000000000000 + + + 0.200000000000000 + + + + + + + + + + + + + Automatic acceleration + + + + + + Variational Bazes Inference acceleration + + + VBI acceleration + + + true + + + + + + + Downsampling and deformation vector interpolation + + + BCPD++ acceleration + + + true + + + + + + + + + + + + + + + + Note: Downsampling automaticallz activates the deformation vecttor interpolation + + + Downsampling + + + true + + + false + + + true + + + + + + Downsampling options. For syntax, please check the documentation + + + Options: + + + + + + + Downsampling options. For syntax, please check the documentation + + + B,5000,0.08 + + + + + + + + + + Convergence + + + true + + + false + + + + + + Tolerance: + + + + + + + 8 + + + 1.000000000000000 + + + 0.000100000000000 + + + 0.000100000000000 + + + + + + + Max number of iterations: + + + + + + + 1 + + + 1000000 + + + 1000 + + + + + + + 1 + + + 1000 + + + 30 + + + + + + + Min number of iterations: + + + + + + + + + + Normalization + + + false + + + + + + Please, refer to the documentation for the explanation of choices + + + Options: + + + + + + + Please, refer to the documentation for the explanation of choices + + + + + + + + + + + + + + + + Postprocessing parameters + + + false + + + + + + Clustering scaling: + + + + + + + 0.000000000000000 + + + 1.000000000000000 + + + + + + + Smoothing iterations: + + + + + + + 10 + + + 1000 + + + + + + + + + + Reset parameters to default + + + + + + + true + + + Generate + + + + + + + + ctkCollapsibleButton + QWidget +
ctkCollapsibleButton.h
+ 1 + + + ctkCollapsibleGroupBox + QGroupBox +
ctkCollapsibleGroupBox.h
+ 1 + + + qMRMLNodeComboBox + QWidget +
qMRMLNodeComboBox.h
+ + + qMRMLWidget + QWidget +
qMRMLWidget.h
+ 1 + + + + + + UI + mrmlSceneChanged(vtkMRMLScene*) + sourceNodeSelectionBox + setMRMLScene(vtkMRMLScene*) + + + 272 + 237 + + + 335 + 57 + + + + + bcpdAdvancedParametersCheckBox + toggled(bool) + bcpdAdvancedControlsGroupBox + setVisible(bool) + + + 114 + 577 + + + 274 + 1065 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelInputMeshLabel + setVisible(bool) + + + 156 + 818 + + + 132 + 847 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelInputMeshLineEdit + setVisible(bool) + + + 156 + 818 + + + 611 + 847 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelNeighboursLabel + setHidden(bool) + + + 156 + 818 + + + 130 + 878 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelNeighboursSpinBox + setHidden(bool) + + + 156 + 818 + + + 611 + 878 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelRadiusLabel + setHidden(bool) + + + 156 + 818 + + + 145 + 910 + + + + + bcpdGeodesicKernelInputMeshCheckBox + toggled(bool) + bcpdGeodesicKernelRadiusDoubleSpinBox + setHidden(bool) + + + 156 + 818 + + + 611 + 910 + + + + + bcpdAccelerationAutomaticPlusPlusCheckBox + toggled(bool) + bcpdDownsamplingCollapsibleGroupBox + setHidden(bool) + + + 134 + 1207 + + + 102 + 1145 + + + + + diff --git a/SlicerBoneMorphing/SlicerBoneMorphing.py b/SlicerBoneMorphing/SlicerBoneMorphing.py new file mode 100644 index 0000000..29673fc --- /dev/null +++ b/SlicerBoneMorphing/SlicerBoneMorphing.py @@ -0,0 +1,22 @@ +from slicer.ScriptedLoadableModule import ScriptedLoadableModule +from src.widget.SlicerBoneMorphingWidget import SlicerBoneMorphingWidget + +class SlicerBoneMorphing(ScriptedLoadableModule): + """Uses ScriptedLoadableModule base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def __init__(self, parent): + ScriptedLoadableModule.__init__(self, parent) + self.parent.title = "Slicer Bone Morphing" + self.parent.categories = ["Morphing"] + self.parent.dependencies = [] + self.parent.contributors = ["Jan Heres (University of West Bohemia), Eva C. Herbst (ETH Zurich), Arthur Porto (Louisiana State University)"] + self.parent.helpText = """ + This module gives a user ability to recreate bone mesh models based on their partial scan. + Please, start with importing your model and checking out the options on the left side afterwards. + """ + self.parent.acknowledgementText = """ + Credits: Jan Heres + I would love to thank my awesome colleagues Eva C. Herbst and Arthur Porto for their priceless contributions to this project! + """ diff --git a/SlicerBoneMorphing/src/CMakeLists.txt b/SlicerBoneMorphing/src/CMakeLists.txt new file mode 100644 index 0000000..1fc302d --- /dev/null +++ b/SlicerBoneMorphing/src/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(logic) +add_subdirectory(widget) diff --git a/SlicerBoneMorphing/src/logic/Constants.py b/SlicerBoneMorphing/src/logic/Constants.py new file mode 100644 index 0000000..5a1073b --- /dev/null +++ b/SlicerBoneMorphing/src/logic/Constants.py @@ -0,0 +1,120 @@ +### RETURN VALUES ### +from src.logic.Enums import BcpdAccelerationMode, BcpdKernelType, BcpdNormalizationOptions, BcpdStandardKernel + +EXIT_OK = 0 +EXIT_FAILURE = 1 + +BCPD_MULTIPLE_VALUES_SEPARATOR = "," +""" + Used when one parameter have multiple values, e.g. -G [string,real,file] +""" + +VALUE_NODE_NOT_SELECTED = 0 # Via Slicer documentation + +### BCPD PARAMETERS### +BCPD_KEY = "bcpd" + +## Tuning ## +BCPD_VALUE_KEY_OMEGA = "-w" +BCPD_VALUE_KEY_LAMBDA = "-l" +BCPD_VALUE_KEY_BETA = "-b" +BCPD_VALUE_KEY_GAMMA = "-g" +BCPD_VALUE_KEY_KAPPA = "-k" + +## Kernel ## +BCPD_VALUE_KEY_KERNEL = "-G" + +## Acceleration ## +BCPD_VALUE_KEY_NYSTORM_G = "-K" +BCPD_VALUE_KEY_NYSTORM_P = "-J" +BCPD_VALUE_KEY_NYSTORM_R = "-r" + +BCPD_VALUE_KEY_KD_TREE = "-p" +BCPD_VALUE_KEY_KD_TREE_SCALE = "-d" +BCPD_VALUE_KEY_KD_TREE_RADIUS = "-e" +BCPD_VALUE_KEY_KD_TREE_THRESHOLD = "-f" + +## Downsampling ## +BCPD_VALUE_KEY_DOWNSAMPLING = "-D" + +## Convergence ## +BCPD_VALUE_KEY_CONVERGENCE_TOLERANCE = "-c" +BCPD_VALUE_KEY_CONVERGENCE_MAX_ITERATIONS = "-n" +BCPD_VALUE_KEY_CONVERGENCE_MIN_ITERATIONS = "-N" + +## Normalization ## +BCPD_VALUE_KEY_NORMALIZATION = "-u" + +## Tuning parameters ## +BCPD_DEFAULT_VALUE_OMEGA = 0.1 +BCPD_DEFAULT_VALUE_LAMBDA = 10 +BCPD_DEFAULT_VALUE_BETA = 10 +BCPD_DEFAULT_VALUE_GAMMA = 0.1 +BCPD_DEFAULT_VALUE_KAPPA = 1000 + +## Kernel functions ## +BCPD_DEFAULT_VALUE_KERNEL_TYPE = BcpdKernelType.STANDARD.value +BCPD_DEFAULT_VALUE_STANDARD_KERNEL_TYPE = BcpdStandardKernel.G0.value +BCPD_DEFAULT_VALUE_TAU = 1 +BCPD_DEFAULT_VALUE_INPUT_MESH_PATH = "" +BCPD_DEFAULT_VALUE_KERNEL_NEIGBOURS = 10 +BCPD_DEFAULT_VALUE_KERNEL_NEIGHBOUR_RADIUS = 10 +BCPD_DEFAULT_VALUE_KERNEL_BETA = 1 +BCPD_DEFAULT_VALUE_KERNEL_K_TILDE = 1 +BCPD_DEFAULT_VALUE_KERNEL_EPSILON = 1 + +## Acceleration mode ## +BCPD_DEFAULT_VALUE_ACCELERATION_MODE = BcpdAccelerationMode.AUTOMATIC.value +BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_G = 140 +BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_J = 600 +BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_R = 1 +BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_SCALE = 7 +BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_RADIUS = 0.3 +BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_SIGMA_THRESHOLD = 0.3 + +## Downsampling ## +BCPD_DEFAULT_VALUE_DOWNSAMPLING_OPTIONS = "B,5000,0.08" + +## Convergence ## +BCPD_DEFAULT_VALUE_CONVERGENCE_TOLERANCE = 0.000001 +BCPD_DEFAULT_VALUE_CONVERGENCE_MAX_ITERATIONS = 1000 +BCPD_DEFAULT_VALUE_CONVERGENCE_MIN_ITERATIONS = 1 + +## Normalization ## +BCPD_DEFAULT_VALUE_NORMALIZATION_OPTIONS = BcpdNormalizationOptions.X.value + +### PREPROCESSING PARAMETERS ### +PREPROCESSING_KEY = "preprocessing" +PREPROCESSING_KEY_DOWNSAMPLING_DISTANCE_THRESHOLD = "ddt" +PREPROCESSING_KEY_NORMALS_ESTIMATION_RADIUS = "ner" +PREPROCESSING_KEY_FPFH_ESTIMATION_RADIUS = "fer" +PREPROCESSING_KEY_MAX_NN_NORMALS = "mnnn" +PREPROCESSING_KEY_MAX_NN_FPFH = "mnf" + +PREPROCESSING_DEFAULT_VALUE_DOWNSAMPLING_DISTANCE_THRESHOLD = 0.1 +PREPROCESSING_DEFAULT_VALUE_RADIUS_NORMAL_SCALE = 4 +PREPROCESSING_DEFAULT_VALUE_RADIUS_FEATURE_SCALE = 10 +PREPROCESSING_DEFAULT_VALUE_MAX_NN_NORMALS = 30 +PREPROCESSING_DEFAULT_VALUE_MAX_NN_FPFH = 100 + +### REGISTRATION PARAMETERS ### +REGISTRATION_KEY_DISTANCE_THRESHOLD = "rdt" +REGISTRATION_KEY_FITNESS_THRESHOLD = "rft" +REGISTRATION_KEY_MAX_ITERATIONS = "rmi" +REGISTRATION_KEY_ICP_DISTANCE_THRESHOLD = "idt" + +REGISTRATION_DEFAULT_VALUE_DISTANCE_THRESHOLD = 1 +REGISTRATION_DEFAULT_VALUE_FITNESS_THRESHOLD = 0.999 +REGISTRATION_DEFAULT_VALUE_MAX_ITERATIONS = 30 +REGISTRATION_DEFAULT_VALUE_ICP_DISTANCE_THRESHOLD = 1 + +### POSTPROCESSING PARAMETERS ### +POSTPROCESSING_KEY = "postprocessing" +POSTPROCESSING_KEY_CLUSTERING_SCALING = "cs" +POSTPROCESSING_KEY_SMOOTHING_ITERATIONS = "sis" + +POSTPROCESSING_DEFAULT_VALUE_CLUSTERING_SCALING = 1 +POSTPROCESSING_DEFAULT_VALUE_SMOOTHING_ITERATIONS = 10 + +### MAX VALUES ### +BCPD_MAX_VALUE_KAPPA = 1000 diff --git a/SlicerBoneMorphing/src/logic/Enums.py b/SlicerBoneMorphing/src/logic/Enums.py new file mode 100644 index 0000000..751d850 --- /dev/null +++ b/SlicerBoneMorphing/src/logic/Enums.py @@ -0,0 +1,21 @@ +from enum import Enum + +class BcpdKernelType(Enum): + STANDARD = 0 + GEODESIC = 1 + +class BcpdStandardKernel(Enum): + G0 = 0 + G1 = 1 + G2 = 2 + G3 = 3 + +class BcpdNormalizationOptions(Enum): + E = 0 + X = 1 + Y = 2 + N = 3 + +class BcpdAccelerationMode(Enum): + AUTOMATIC = 0 + MANUAL = 1 diff --git a/SlicerBoneMorphing/src/logic/SlicerBoneMorphingLogic.py b/SlicerBoneMorphing/src/logic/SlicerBoneMorphingLogic.py new file mode 100644 index 0000000..8162bf1 --- /dev/null +++ b/SlicerBoneMorphing/src/logic/SlicerBoneMorphingLogic.py @@ -0,0 +1,460 @@ +from sys import platform +from typing import Tuple +from slicer.ScriptedLoadableModule import ScriptedLoadableModuleLogic +import slicer +from slicer import vtkMRMLModelNode + +try: + import open3d as o3d +except ModuleNotFoundError: + print("Module Open3D is not installed. Installing...") + slicer.util.pip_install('open3d') + +import open3d as o3d +import numpy as np +import os +import glob +import vtk +from vtk.util.numpy_support import vtk_to_numpy +import subprocess +import tempfile + +from .Constants import * + +BCPD_EXEC = os.path.dirname( + os.path.abspath(__file__)) + "/../../Resources/BCPD/exec/" + +# NOTE: Needs relative path to the main module script +if platform == "linux" or platform == "linux2": + BCPD_EXEC += "bcpd_linux_x86_64" +elif platform == "darwin": + # Slicer is running through Rosetta, so x86 version needs to be used for now + BCPD_EXEC += "bcpd_macos_x86_64" +elif platform == "win32": + BCPD_EXEC += "bcpd_win32.exe" + +deformed = None + + +class SlicerBoneMorphingLogic(ScriptedLoadableModuleLogic): + """This class should implement all the actual + computation done by your module. The interface + should be such that other python code can import + this class and make use of the functionality without + requiring an instance of the Widget. + Uses ScriptedLoadableModuleLogic base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def __init__(self, parent): + ScriptedLoadableModuleLogic.__init__(self, parent) + + def generate_model( + self, source_model: vtkMRMLModelNode, + target_model: vtkMRMLModelNode, + parameters) -> Tuple[int, vtk.vtkPolyData, vtk.vtkPolyData]: + """ + Generates new model based on the BCPD algorithm fit between source and target models. + + Parameters + ---------- + vtkMRMLModelNode source_model - source (partial) model to be fit-generated + vtkMRMLModelNode target_model - model to fit the partial source by. + string parameters - parameters for the preprocessing, BCPD and postprocessing + + Returns + ------- + Tuple[int status, vtk.vtkPolyData generatedPolydata, vtk.vtkPolyData mergedPolydata]: + - status: EXIT_OK or EXIT_FAILURE + - generatedPolydata: Generated model by the BCPD + - mergedPolydata: generatedPolydata that had been merged with the targetModel + """ + + if (source_model == VALUE_NODE_NOT_SELECTED + or target_model == VALUE_NODE_NOT_SELECTED): + print("Input or foundation model(s) were not selected") + return EXIT_FAILURE, None, None + + source_mesh = self.convert_model_to_mesh(source_model) + target_mesh = self.convert_model_to_mesh(target_model) + + err, result_icp = self.preprocess_model(source_mesh, target_mesh, + parameters[PREPROCESSING_KEY]) + if err == EXIT_FAILURE: + print("Cannot continue with generating. Aborting...") + return EXIT_FAILURE, None, None + + # BCPD stage + deformed = self.deformable_registration( + source_mesh.transform(result_icp.transformation), target_mesh, + parameters[BCPD_KEY]) + if (deformed == None): + return EXIT_FAILURE, None, None + + generated_polydata, merged_polydata = self.postprocess_data( + deformed, target_mesh, parameters[POSTPROCESSING_KEY]) + + return EXIT_OK, generated_polydata, merged_polydata + + def convert_mesh_to_vtk_polydata( + self, mesh: o3d.geometry.TriangleMesh) -> vtk.vtkPolyData: + """ + Convert o3d.geometry.TriangleMesh to vtkPolyData + + Parameters + ---------- + o3d.geometry.TriangleMesh mesh - mesh to be converted + + Returns + ------- + vtk.vtkPolyData representation of the input mesh + """ + vertices = np.asarray(mesh.vertices) + triangles = np.asarray(mesh.triangles) + + points = vtk.vtkPoints() + + for v in vertices: + points.InsertNextPoint(v) + + polys = vtk.vtkCellArray() + + for t in triangles: + polys.InsertNextCell(3, t) + + vtk_poly_data = vtk.vtkPolyData() + vtk_poly_data.SetPoints(points) + vtk_poly_data.SetPolys(polys) + + return vtk_poly_data + + def convert_model_to_mesh(self, model: vtkMRMLModelNode): + """ + Convert vtkMRMLModelNode to open3d.geometry.TriangleMesh + + Parameters + ---------- + slicer.vtkMRMLNode model - model to be converted + + Returns + ------- + open3d.geometry.TriangleMesh + """ + vtk_poly_data = model.GetPolyData() + + # Get vertices from vtk_polydata + numpy_vertices = vtk_to_numpy(vtk_poly_data.GetPoints().GetData()) + + # Get normals if present + numpy_normals = None + model_normals = vtk_poly_data.GetPointData().GetNormals() + if (model_normals is not None): + numpy_normals = vtk_to_numpy( + vtk_poly_data.GetPointData().GetNormals()) + + # Get indices (triangles), this would be a (n, 3) shape numpy array where n is the number of triangles. + # If the vtkPolyData does not represent a valid triangle mesh, the indices may not form valid triangles. + numpy_indices = vtk_to_numpy( + vtk_poly_data.GetPolys().GetData()).reshape(-1, 4)[:, 1:] + + # convert numpy array to open3d TriangleMesh + open3d_mesh = o3d.geometry.TriangleMesh( + o3d.utility.Vector3dVector(numpy_vertices), + o3d.utility.Vector3iVector(numpy_indices)) + + if numpy_normals is not None: + open3d_mesh.vertex_normals = o3d.utility.Vector3dVector( + numpy_normals) + + return open3d_mesh + + def convert_to_point_cloud( + self, mesh: o3d.geometry.TriangleMesh) -> o3d.geometry.PointCloud: + """ + Convert o3d.geometry.TriangleMesh to o3d.geometry.PointCloud + + Parameters + ---------- + o3d.geometry.TriangleMesh mesh - mesh to be converted + + Returns + ------- + o3d.geometry.PointCloud + """ + + # mesh_center = mesh.get_center() + # mesh.translate(-mesh_center, relative=False) # Not needed for Slicer + pcd = o3d.geometry.PointCloud() + pcd.points = mesh.vertices + pcd.colors = mesh.vertex_colors + pcd.normals = mesh.vertex_normals + + return pcd + + def preprocess_model( + self, source_mesh: o3d.geometry.TriangleMesh, + target_mesh: o3d.geometry.TriangleMesh, parameters: dict + ) -> Tuple[int, o3d.pipelines.registration.RegistrationResult]: + """ + Preprocess model before advancing into the generation (BCPD) stage. This method converts the input models into respective point clouds, then performs RANSAC best fit transformation from source to target and returns the result + + Parameters + ---------- + o3d.geometry.TriangleMesh source_mesh: Source MRML model + o3d.geometry.TriangleMesh target_mesh: Target MRML model + + Returns + ------- + Tuple[int, o3d.pipelines.registration.RegistrationResult] - if Tuple[0] equals EXIT_OK, then Tuple[1] will carry the registration result + """ + + source_pcd = self.convert_to_point_cloud(source_mesh) + target_pcd = self.convert_to_point_cloud(target_mesh) + + source_pcd_downsampled, source_pcd_fpfh = self.preprocess_point_cloud( + source_pcd, + parameters[PREPROCESSING_KEY_DOWNSAMPLING_DISTANCE_THRESHOLD], + parameters[PREPROCESSING_KEY_NORMALS_ESTIMATION_RADIUS], + parameters[PREPROCESSING_KEY_FPFH_ESTIMATION_RADIUS], + parameters[PREPROCESSING_KEY_MAX_NN_NORMALS], + parameters[PREPROCESSING_KEY_MAX_NN_FPFH]) + + target_pcd_downsampled, target_pcd_fpfh = self.preprocess_point_cloud( + target_pcd, + parameters[PREPROCESSING_KEY_DOWNSAMPLING_DISTANCE_THRESHOLD], + parameters[PREPROCESSING_KEY_NORMALS_ESTIMATION_RADIUS], + parameters[PREPROCESSING_KEY_FPFH_ESTIMATION_RADIUS], + parameters[PREPROCESSING_KEY_MAX_NN_NORMALS], + parameters[PREPROCESSING_KEY_MAX_NN_FPFH]) + + try: + result_ransac = self.ransac_pcd_registration( + source_pcd_downsampled, target_pcd_downsampled, + source_pcd_fpfh, target_pcd_fpfh, + parameters[REGISTRATION_KEY_DISTANCE_THRESHOLD], + parameters[REGISTRATION_KEY_FITNESS_THRESHOLD], + parameters[REGISTRATION_KEY_MAX_ITERATIONS]) + except RuntimeError: + print( + "No ideal fit was found using the RANSAC algorithm. Please, try adjusting the parameters" + ) + return EXIT_FAILURE, None + + result_icp = o3d.pipelines.registration.registration_icp( + source_pcd_downsampled, target_pcd_downsampled, + parameters[REGISTRATION_KEY_ICP_DISTANCE_THRESHOLD], + result_ransac.transformation, + o3d.pipelines.registration.TransformationEstimationPointToPlane()) + + return EXIT_OK, result_icp + + def preprocess_point_cloud( + self, pcd: o3d.geometry.PointCloud, + downsampling_distance_threshold: float, + normals_estimation_radius: float, fpfh_estimation_radius: float, + max_nn_normals: int, max_nn_fpfh: int + ) -> Tuple[o3d.geometry.PointCloud, o3d.pipelines.registration.Feature]: + ''' + Perform downsampling of a mesh, normal estimation and computing FPFH feature of the point cloud. + + Parameters + ---------- + o3d.geometry.PointCloud pcd: Source point cloud + float downsampling_distance_threshold: Distance threshold for downsampling + float normals_estimation_radius: Radius for estimating normals + float fpfh_estimation_radius: Radius for the FPFH computation + int max_nn_normals: Maximum number of neighbours considered for normals estimation + int max_nn_fpfh: Maximum number of neighbours considered for the FPFH calculation + + Returns + ------- + Tuple: [Downsampled PCD: open3d.geometry.PointCloud, + FPFH: open3d.pipelines.registration.Feature] + ''' + + pcd_downsampled: o3d.geometry.PointCloud = pcd.voxel_down_sample( + downsampling_distance_threshold) + + pcd_downsampled.estimate_normals( + o3d.geometry.KDTreeSearchParamHybrid( + radius=normals_estimation_radius, max_nn=max_nn_normals)) + + pcd_fpfh = o3d.pipelines.registration.compute_fpfh_feature( + pcd_downsampled, + o3d.geometry.KDTreeSearchParamHybrid(radius=fpfh_estimation_radius, + max_nn=max_nn_fpfh)) + + return pcd_downsampled, pcd_fpfh + + def ransac_pcd_registration( + self, + source_pcd_down: o3d.geometry.PointCloud, + target_pcd_down: o3d.geometry.PointCloud, + source_fpfh: o3d.pipelines.registration.Feature, + target_fpfh: o3d.pipelines.registration.Feature, + distance_threshold: float, + fitness_threshold: float, + max_iterations: int, + ) -> o3d.pipelines.registration.RegistrationResult: + ''' Perform a registration of nearest neighbours using the RANSAC algorithm. + + Parameters + ---------- + o3d.geometry.PointCloud source_pcd_down: Downsampled SOURCE point cloud + o3d.geometry.PointCloud target_pcd_down: Downsampled TARGET point cloud + o3d.pipelines.registration.Feature source_fpfh: Source PCD Fast-Point-Feature-Histogram + o3d.pipelines.registration.Feature target_fpfh: Target PCD Fast-Point-Feature-Histogram + float distance_threshold: Threshold in which a near point is considered a neighbour + float fitness_threshold: Minimal value for iterations until it is reached + max_iterations: Maximum number of iterations of the RANSAC algorithm + + Returns + ------- + Best PCD fit: open3d.pipelines.registration.RegistrationResult + + ''' + fitness = 0 + count = 0 + best_result = None + fitness_max = 1 + + while (fitness < fitness_threshold and fitness < fitness_max + and count < max_iterations): + result = o3d.pipelines.registration.registration_ransac_based_on_feature_matching( + source_pcd_down, + target_pcd_down, + source_fpfh, + target_fpfh, + True, + distance_threshold, + o3d.pipelines.registration. + TransformationEstimationPointToPoint(True), + 3, + [ + o3d.pipelines.registration. + CorrespondenceCheckerBasedOnEdgeLength(0.9), + o3d.pipelines.registration. + CorrespondenceCheckerBasedOnDistance(distance_threshold) + ], + # NOTE: Just for earlier termination, but still needs the outer loop for proper convergence + o3d.pipelines.registration.RANSACConvergenceCriteria( + 100000, fitness_threshold)) + if result.fitness > fitness and result.fitness < 1: + fitness = result.fitness + best_result = result + count += 1 + return best_result + + def deformable_registration(self, source_pcd: o3d.geometry.PointCloud, + target_pcd: o3d.geometry.PointCloud, + bcpd_parameters) -> o3d.geometry.TriangleMesh: + """ + Perform a BCPD deformable registration from the source point cloud to the target point cloud + + Parameters + ---------- + o3d.geometry.PointCloud sourcePcd: source point cloud + o3d.geometry.PointCloud targetPcd: target point cloud + string bcpdParameters: parameters for the BCPD algorithm + + Returns + ------- + o3d.geometry.TriangleMesh representing the new deformed mesh + """ + source_array = np.asarray(source_pcd.vertices, dtype=np.float32) + target_array = np.asarray(target_pcd.vertices, dtype=np.float32) + + target_path = tempfile.gettempdir( + ) + '/slicer_bone_morphing_target.txt' + source_path = tempfile.gettempdir( + ) + '/slicer_bone_morphing_source.txt' + output_path = tempfile.gettempdir() + '/output_' + + np.savetxt(target_path, target_array, delimiter=',') + np.savetxt(source_path, source_array, delimiter=',') + + cmd = f'{BCPD_EXEC} -h -x {target_path} -y {source_path}' + + for key in bcpd_parameters.keys(): + cmd += f' {key}{bcpd_parameters[key]}' + + cmd += f' -o {output_path}' + print("BCPD: " + cmd) + + subprocess.run(cmd, + shell=True, + check=True, + text=True, + capture_output=True) + + try: + bcpdResult = np.loadtxt(output_path + "y.interpolated.txt") + except FileNotFoundError: + print( + "No results generated by BCPD. Refer to the output in the console." + ) + return None + + for fl in glob.glob(output_path + "*.txt"): + os.remove(fl) + os.remove(target_path) + os.remove(source_path) + + deformed = o3d.geometry.TriangleMesh() + deformed.vertices = o3d.utility.Vector3dVector(np.asarray(bcpdResult)) + deformed.triangles = source_pcd.triangles + + return deformed + + def postprocess_data( + self, deformed: o3d.geometry.TriangleMesh, + target_mesh: o3d.geometry.TriangleMesh, + parameters: dict) -> Tuple[vtk.vtkPolyData, vtk.vtkPolyData]: + """ + Perform a postprocessing (smoothing and filtering) of the result of the BCPD algorithm + + Parameters + ---------- + o3d.geometry.TriangleMesh deformed: BCPD result mesh + o3d.geometry.TriangleMesh target_mesh: Target model mesh + dict parameters: parameters dict for the postprocessing stage + + Returns + ------- + Tuple[vtk.vtkPolyData, vtk.vtkPolyData]: Tuple[0] will represent the postprocessed BCPD mesh, Tuple[1] will represent the COMBINED mesh (i.e. BCPD result merged with the target mesh) + + """ + deformed.compute_vertex_normals() + target_mesh.compute_vertex_normals() + + # Combine meshes (alternative - to crop the first before merging) + combined = deformed + target_mesh + combined.compute_vertex_normals() + + # Simplify mesh (smoothing and filtering) + deformed_smp = deformed.simplify_vertex_clustering( + parameters[POSTPROCESSING_KEY_CLUSTERING_SCALING], + contraction=o3d.geometry.SimplificationContraction.Average) + deformed_smp = deformed_smp.filter_smooth_simple( + number_of_iterations=parameters[ + POSTPROCESSING_KEY_SMOOTHING_ITERATIONS]) + deformed_smp = deformed_smp.filter_smooth_taubin( + number_of_iterations=parameters[ + POSTPROCESSING_KEY_SMOOTHING_ITERATIONS]) + + mesh_smp = combined.simplify_vertex_clustering( + parameters[POSTPROCESSING_KEY_CLUSTERING_SCALING], + contraction=o3d.geometry.SimplificationContraction.Average) + mesh_smp = mesh_smp.filter_smooth_simple( + number_of_iterations=parameters[ + POSTPROCESSING_KEY_SMOOTHING_ITERATIONS]) + mesh_smp = mesh_smp.filter_smooth_taubin( + number_of_iterations=parameters[ + POSTPROCESSING_KEY_SMOOTHING_ITERATIONS]) + + # mesh_smp.compute_vertex_normals() # NOTE: Is this needed? + + generated_polydata = self.convert_mesh_to_vtk_polydata(deformed_smp) + merged_polydata = self.convert_mesh_to_vtk_polydata(mesh_smp) + + return generated_polydata, merged_polydata diff --git a/SlicerBoneMorphing/src/widget/SlicerBoneMorphingWidget.py b/SlicerBoneMorphing/src/widget/SlicerBoneMorphingWidget.py new file mode 100644 index 0000000..9d5c323 --- /dev/null +++ b/SlicerBoneMorphing/src/widget/SlicerBoneMorphingWidget.py @@ -0,0 +1,340 @@ +from src.logic.Constants import * +import ctk +import slicer +from slicer.ScriptedLoadableModule import ScriptedLoadableModuleWidget +from src.logic.SlicerBoneMorphingLogic import SlicerBoneMorphingLogic + +from src.logic.Enums import BcpdAccelerationMode, BcpdKernelType, BcpdNormalizationOptions, BcpdStandardKernel +from qt import QComboBox + +import os + + +class SlicerBoneMorphingWidget(ScriptedLoadableModuleWidget): + """Uses ScriptedLoadableModuleWidget base class, available at: + https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py + """ + + def __init__(self, parent): + """Called when the application opens the module the first time and the widget is initialized.""" + ScriptedLoadableModuleWidget.__init__(self, parent) + self.bcpd_options = {} + + def setup(self): + """Called when the application opens the module the first time and the widget is initialized.""" + ScriptedLoadableModuleWidget.setup(self) + + # Load widget from .ui file (created by Qt Designer) + self.uiWidget = slicer.util.loadUI( + self.resourcePath("UI/SlicerBoneMorphing.ui")) + self.layout.addWidget(self.uiWidget) + self.ui = slicer.util.childWidgetVariables(self.uiWidget) + self.logic = SlicerBoneMorphingLogic(self) + + self.setup_ui() + self.reset_parameters_to_default() + + def setup_ui(self): + """ + Method that sets up all UI elements and their dependencies + """ + + self.ui.sourceNodeSelectionBox.setMRMLScene(slicer.mrmlScene) + self.ui.targetNodeSelectionBox.setMRMLScene(slicer.mrmlScene) + + self.ui.bcpdAdvancedControlsGroupBox.setVisible(False) + + self.setup_combo_box(self.ui.bcpdKernelTypeComboBox, BcpdKernelType, + self.show_kernel_type) + self.ui.bcpdGeodesicKernelGroupBox.setVisible(False) + + self.setup_combo_box(self.ui.bcpdStandardKernelComboBox, + BcpdStandardKernel, None) + + self.setup_combo_box(self.ui.bcpdAccelerationModeComboBox, + BcpdAccelerationMode, self.show_acceleration_mode) + self.ui.bcpdAccelerationManualGroupBox.setVisible(False) + + self.ui.bcpdGeodesicKernelInputMeshLineEdit.setVisible(False) + self.ui.bcpdGeodesicKernelInputMeshLabel.setVisible(False) + + self.setup_combo_box(self.ui.bcpdNormalizationComboBox, + BcpdNormalizationOptions, None) + + self.ui.bcpdDownsamplingCollapsibleGroupBox.visible = False + + self.ui.bcpdResetParametersPushButton.clicked.connect( + self.reset_parameters_to_default) + self.ui.generateModelButton.clicked.connect(self.generate_model) + + def reset_parameters_to_default(self): + ## Preprocessing parameters ## + self.ui.preprocessingDownsamplingDistanceThresholdDoubleSpinBox.value = PREPROCESSING_DEFAULT_VALUE_DOWNSAMPLING_DISTANCE_THRESHOLD + self.ui.preprocessingNormalsEstimationRadiusDoubleSpinBox.value = PREPROCESSING_DEFAULT_VALUE_RADIUS_NORMAL_SCALE + self.ui.preprocessingNormalsEstimationMaxNeighboursSpinBox.value = PREPROCESSING_DEFAULT_VALUE_MAX_NN_NORMALS + self.ui.preprocessingFpfhRadiusDoubleSpinBox.value = PREPROCESSING_DEFAULT_VALUE_RADIUS_FEATURE_SCALE + self.ui.preprocessingFpfhMaxNeighboursSpinBox.value = PREPROCESSING_DEFAULT_VALUE_MAX_NN_FPFH + + ## Registration parameters ## + self.ui.registrationMaxIterationsSpinBox.value = REGISTRATION_DEFAULT_VALUE_MAX_ITERATIONS + self.ui.registrationDistanceThresholdDoubleSpinBox.value = REGISTRATION_DEFAULT_VALUE_DISTANCE_THRESHOLD + self.ui.registrationFitnessThresholdDoubleSpinBox.value = REGISTRATION_DEFAULT_VALUE_FITNESS_THRESHOLD + self.ui.registrationIcpDistanceThresholdDoubleSpinBox.value = REGISTRATION_DEFAULT_VALUE_ICP_DISTANCE_THRESHOLD + + ## Tuning parameters ## + self.ui.bcpdOmegaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_OMEGA + self.ui.bcpdLambdaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_LAMBDA + self.ui.bcpdBetaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_BETA + self.ui.bcpdGammaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_GAMMA + self.ui.bcpdKappaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_KAPPA + + ## Kernel parameters ## + self.ui.bcpdKernelTypeComboBox.setCurrentIndex( + BCPD_DEFAULT_VALUE_KERNEL_TYPE) + self.ui.bcpdStandardKernelComboBox.setCurrentIndex( + BCPD_DEFAULT_VALUE_STANDARD_KERNEL_TYPE) + self.ui.bcpdGeodesicKernelTauDoubleSpinBox.value = BCPD_DEFAULT_VALUE_TAU + self.ui.bcpdGeodesicKernelInputMeshCheckBox.checked = False + self.ui.bcpdGeodesicKernelInputMeshLineEdit.text = BCPD_DEFAULT_VALUE_INPUT_MESH_PATH + self.ui.bcpdGeodesicKernelNeighboursSpinBox.value = BCPD_DEFAULT_VALUE_KERNEL_NEIGBOURS + self.ui.bcpdGeodesicKernelRadiusDoubleSpinBox.value = BCPD_DEFAULT_VALUE_KERNEL_NEIGHBOUR_RADIUS + self.ui.bcpdGeodesicKernelBetaDoubleSpinBox.value = BCPD_DEFAULT_VALUE_KERNEL_BETA + self.ui.bcpdGeodesicKernelKTildeDoubleSpinBox.value = BCPD_DEFAULT_VALUE_KERNEL_K_TILDE + self.ui.bcpdGeodesicKernelEpsilonDoubleSpinBox.value = BCPD_DEFAULT_VALUE_KERNEL_EPSILON + + ## Acceleration parameters ## + self.ui.bcpdAccelerationModeComboBox.setCurrentIndex( + BCPD_DEFAULT_VALUE_ACCELERATION_MODE) + self.ui.bcpdAccelerationAutomaticVbiCheckBox.checked = True + self.ui.bcpdAccelerationAutomaticPlusPlusCheckBox.checked = True + self.ui.bcpdAccelerationManualNystormGSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_G + self.ui.bcpdAccelerationManualNystormJSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_J + self.ui.bcpdAccelerationManualNystormRSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_NYSTORM_SAMPLES_R + self.ui.bcpdAccelerationManualKdTreeScaleDoubleSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_SCALE + self.ui.bcpdAccelerationManualKdTreeRadiusDoubleSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_RADIUS + self.ui.bcpdAccelerationManualKdTreeThresholdDoubleSpinBox.value = BCPD_DEFAULT_VALUE_ACCELERATION_KD_TREE_SIGMA_THRESHOLD + + ## Downsampling options ## + self.ui.bcpdDownsamplingLineEdit.text = BCPD_DEFAULT_VALUE_DOWNSAMPLING_OPTIONS + + ## Convergence options ## + self.ui.bcpdConvergenceToleranceDoubleSpinBox.value = BCPD_DEFAULT_VALUE_CONVERGENCE_TOLERANCE + self.ui.bcpdConvergenceMaxIterationsSpinBox.value = BCPD_DEFAULT_VALUE_CONVERGENCE_MAX_ITERATIONS + self.ui.bcpdConvergenceMinIterationsSpinBox.value = BCPD_DEFAULT_VALUE_CONVERGENCE_MIN_ITERATIONS + + ## Normalization options ## + self.ui.bcpdNormalizationComboBox.setCurrentIndex( + BCPD_DEFAULT_VALUE_NORMALIZATION_OPTIONS) + + self.ui.postprocessingClusteringScalingDoubleSpinBox.value = POSTPROCESSING_DEFAULT_VALUE_CLUSTERING_SCALING + self.ui.processingSmoothingIterationsSpinBox.value = POSTPROCESSING_DEFAULT_VALUE_SMOOTHING_ITERATIONS + + def setup_combo_box(self, combo_box: QComboBox, enum, onSelectionChanged): + """ + Method for setting up combo box and its possible values + """ + for mode in list(enum): + combo_box.addItem(mode.name, mode) + if onSelectionChanged is not None: + combo_box.currentIndexChanged.connect(onSelectionChanged) + combo_box.setCurrentIndex(0) + + def show_kernel_type(self, current_index): + """ + Kernel type callback + """ + if current_index == BcpdKernelType.STANDARD.value: + show_standard_setting = True + show_geodesic_settings = False + else: + show_standard_setting = False + show_geodesic_settings = True + + self.ui.bcpdStandardKernelGroupBox.setVisible(show_standard_setting) + self.ui.bcpdGeodesicKernelGroupBox.setVisible(show_geodesic_settings) + + def show_acceleration_mode(self, currentIndex): + """ + Acceleration mode combo box callback + """ + if currentIndex == BcpdAccelerationMode.AUTOMATIC.value: + show_automatic = True + show_manual = False + else: + show_automatic = False + show_manual = True + + self.ui.bcpdAccelerationAutomaticGroupBox.setVisible(show_automatic) + self.ui.bcpdAccelerationManualGroupBox.setVisible(show_manual) + + def parse_parameters(self) -> dict: + """ + Parsing parameters from the UI user option elements + """ + params = {} + params[PREPROCESSING_KEY] = self.parse_parameters_preprocessing() + params[BCPD_KEY] = self.parse_parameters_bcpd() + params[POSTPROCESSING_KEY] = self.parse_parameters_postprocessing() + + return params + + def parse_parameters_preprocessing(self) -> dict: + params = {} + + ## Preprocessing + params[ + PREPROCESSING_KEY_DOWNSAMPLING_DISTANCE_THRESHOLD] = self.ui.preprocessingDownsamplingDistanceThresholdDoubleSpinBox.value + params[ + PREPROCESSING_KEY_NORMALS_ESTIMATION_RADIUS] = self.ui.preprocessingNormalsEstimationRadiusDoubleSpinBox.value + params[ + PREPROCESSING_KEY_MAX_NN_NORMALS] = self.ui.preprocessingNormalsEstimationMaxNeighboursSpinBox.value + params[ + PREPROCESSING_KEY_FPFH_ESTIMATION_RADIUS] = self.ui.preprocessingFpfhRadiusDoubleSpinBox.value + params[ + PREPROCESSING_KEY_MAX_NN_FPFH] = self.ui.preprocessingFpfhMaxNeighboursSpinBox.value + + ## Registration + params[ + REGISTRATION_KEY_MAX_ITERATIONS] = self.ui.registrationMaxIterationsSpinBox.value + params[ + REGISTRATION_KEY_DISTANCE_THRESHOLD] = self.ui.registrationDistanceThresholdDoubleSpinBox.value + params[ + REGISTRATION_KEY_FITNESS_THRESHOLD] = self.ui.registrationFitnessThresholdDoubleSpinBox.value + params[ + REGISTRATION_KEY_ICP_DISTANCE_THRESHOLD] = self.ui.registrationIcpDistanceThresholdDoubleSpinBox.value + + return params + + def parse_parameters_bcpd(self) -> dict: + """ + Parsing parameters from the UI for the BCPD stage + """ + params = {} + + ## Tuning parameters ## + params[BCPD_VALUE_KEY_OMEGA] = self.ui.bcpdOmegaDoubleSpinBox.value + params[BCPD_VALUE_KEY_LAMBDA] = self.ui.bcpdLambdaDoubleSpinBox.value + params[BCPD_VALUE_KEY_BETA] = self.ui.bcpdBetaDoubleSpinBox.value + params[BCPD_VALUE_KEY_GAMMA] = self.ui.bcpdGammaDoubleSpinBox.value + + kappa = self.ui.bcpdKappaDoubleSpinBox.value + if kappa < BCPD_MAX_VALUE_KAPPA: # Setting it to max behaves like "infinity" + params[BCPD_VALUE_KEY_KAPPA] = kappa + + # if self.ui.bcpdAdvancedParametersCheckBox.checked == True: + self.parse_advanced_parameters(params) + + return params + + def parse_advanced_parameters(self, params: dict) -> None: + """ + Parsing parameters under the "Advanced" section + """ + ## Kernel settings ## + kernel_params = "" + kernel_type = self.ui.bcpdKernelTypeComboBox.currentIndex + if (kernel_type == BcpdKernelType.STANDARD.value): + selected_kernel = self.ui.bcpdStandardKernelComboBox.currentIndex + # Default kernel is Gauss, which does not need to be specified + if selected_kernel != BCPD_DEFAULT_VALUE_STANDARD_KERNEL_TYPE: + kernel_params += str(selected_kernel) + else: # Geodesic Kernel + kernel_params += "geodesic" + BCPD_MULTIPLE_VALUES_SEPARATOR + + tau = self.ui.bcpdGeodesicKernelTauDoubleSpinBox.value + kernel_params += str(tau) + BCPD_MULTIPLE_VALUES_SEPARATOR + + if self.ui.bcpdGeodesicKernelInputMeshCheckBox.checked == True: + input_mesh_path = self.ui.bcpdGeodesicKernelInputMeshLineEdit.text + if not os.path.exists(input_mesh_path): + print("File '" + input_mesh_path + + "' does not exist. Cancelling process...") + return + kernel_params += input_mesh_path + else: + kernel_params += str(self.ui.bcpdGeodesicKernelNeighboursSpinBox.value) + \ + BCPD_MULTIPLE_VALUES_SEPARATOR + kernel_params += str( + self.ui.bcpdGeodesicKernelRadiusDoubleSpinBox.value) + + if kernel_params != "": + params[BCPD_VALUE_KEY_KERNEL] = kernel_params + + ## Acceleration settings ## + if self.ui.bcpdAccelerationModeComboBox.currentIndex == BcpdAccelerationMode.AUTOMATIC.value: + if self.ui.bcpdAccelerationAutomaticVbiCheckBox.checked == True: + params[BCPD_VALUE_KEY_NYSTORM_G] = 70 + params[BCPD_VALUE_KEY_NYSTORM_P] = 300 + # Option switch without a value + params[BCPD_VALUE_KEY_KD_TREE] = "" + params[BCPD_VALUE_KEY_KD_TREE_SCALE] = 7 + params[BCPD_VALUE_KEY_KD_TREE_RADIUS] = 0.15 + if self.ui.bcpdAccelerationAutomaticPlusPlusCheckBox.checked == True: + params[BCPD_VALUE_KEY_DOWNSAMPLING] = "B,10000,0.08" + else: # Manual acceleration + if self.ui.bcpdAccelerationManualNystormGroupBox.checked == True: + params[ + BCPD_VALUE_KEY_NYSTORM_G] = self.ui.bcpdAccelerationManualNystormGSpinBox.value + params[ + BCPD_VALUE_KEY_NYSTORM_P] = self.ui.bcpdAccelerationManualNystormJSpinBox.value + params[ + BCPD_VALUE_KEY_NYSTORM_R] = self.ui.bcpdAccelerationManualNystormRSpinBox.value + if self.ui.bcpdAccelerationManualKdTreeGroupBox.checked == True: + # Option switch without a value + params[BCPD_VALUE_KEY_KD_TREE] = "" + params[ + BCPD_VALUE_KEY_KD_TREE_SCALE] = self.ui.bcpdAccelerationManualKdTreeScaleDoubleSpinBox.value + params[ + BCPD_VALUE_KEY_KD_TREE_RADIUS] = self.ui.bcpdAccelerationManualKdTreeRadiusDoubleSpinBox.value + params[ + BCPD_VALUE_KEY_KD_TREE_THRESHOLD] = self.ui.bcpdAccelerationManualKdTreeThresholdDoubleSpinBox.value + + ## Downsampling settings ## + + if params.get(BCPD_VALUE_KEY_DOWNSAMPLING) is None: + params[ + BCPD_VALUE_KEY_DOWNSAMPLING] = self.ui.bcpdDownsamplingLineEdit.text + + ## Convergence options ## + params[ + BCPD_VALUE_KEY_CONVERGENCE_TOLERANCE] = self.ui.bcpdConvergenceToleranceDoubleSpinBox.value + params[ + BCPD_VALUE_KEY_CONVERGENCE_MIN_ITERATIONS] = self.ui.bcpdConvergenceMinIterationsSpinBox.value + params[ + BCPD_VALUE_KEY_CONVERGENCE_MAX_ITERATIONS] = self.ui.bcpdConvergenceMaxIterationsSpinBox.value + + ## Normalization options ## + params[ + BCPD_VALUE_KEY_NORMALIZATION] = self.ui.bcpdNormalizationComboBox.currentText.lower( + ) + + def parse_parameters_postprocessing(self) -> dict: + params = {} + + params[ + POSTPROCESSING_KEY_CLUSTERING_SCALING] = self.ui.postprocessingClusteringScalingDoubleSpinBox.value + params[ + POSTPROCESSING_KEY_SMOOTHING_ITERATIONS] = self.ui.processingSmoothingIterationsSpinBox.value + + return params + + def generate_model(self) -> None: + """ + Generate button callback. Calls the Logic's generate_model method and adds the results into the scene + """ + params = self.parse_parameters() + + err, generated_polydata, merged_polydata = self.logic.generate_model( + self.ui.sourceNodeSelectionBox.currentNode(), + self.ui.targetNodeSelectionBox.currentNode(), params) + + if (err == EXIT_OK): + model_node = slicer.mrmlScene.AddNewNodeByClass( + 'vtkMRMLModelNode', 'BCPD generated') + model_node.SetAndObservePolyData(generated_polydata) + model_node.CreateDefaultDisplayNodes() + + model_node = slicer.mrmlScene.AddNewNodeByClass( + 'vtkMRMLModelNode', 'BCPD merged') + model_node.SetAndObservePolyData(merged_polydata) + model_node.CreateDefaultDisplayNodes() diff --git a/docs/.DS_Store b/docs/.DS_Store new file mode 100644 index 0000000..53ff113 Binary files /dev/null and b/docs/.DS_Store differ diff --git a/docs/assets/results.png b/docs/assets/results.png new file mode 100644 index 0000000..85ca712 Binary files /dev/null and b/docs/assets/results.png differ diff --git a/docs/assets/ui.png b/docs/assets/ui.png new file mode 100644 index 0000000..6fc0b8e Binary files /dev/null and b/docs/assets/ui.png differ diff --git a/metadata.s4ext b/metadata.s4ext new file mode 100644 index 0000000..ec0cf29 --- /dev/null +++ b/metadata.s4ext @@ -0,0 +1,5 @@ +# This is a comment +metadataname This is the value associated with 'metadataname' + +# This is an other comment +anothermetadata This is the value associated with 'anothermetadata'