Skip to content

Commit

Permalink
ENH: Hierarchical approach to 3d registration
Browse files Browse the repository at this point in the history
  • Loading branch information
NicerNewerCar committed Jun 5, 2024
1 parent 7a2dd0f commit b3a8fc1
Show file tree
Hide file tree
Showing 12 changed files with 962 additions and 62 deletions.
66 changes: 57 additions & 9 deletions AutoscoperM/AutoscoperM.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,13 @@
import shutil
import time
import zipfile
from itertools import product
from typing import Optional, Union

import qt
import slicer
import vtk
from numpy.typing import NDArray
from slicer.ScriptedLoadableModule import (
ScriptedLoadableModule,
ScriptedLoadableModuleLogic,
Expand Down Expand Up @@ -1071,15 +1073,12 @@ def progressCallback(x):
tfm.SetElement(1, 3, origin[1])
tfm.SetElement(2, 3, origin[2])

IO.createTRAFile(
volName=segmentName,
trialName=None,
outputDir=outputDir,
trackingsubDir=trackingSubDir,
volSize=segmentVolumeSize,
Origin2DicomTransformFile=origin2DicomTransformFile,
transform=tfm,
)
if origin2DicomTransformFile is not None:
tfm = self.applyOrigin2DicomTransform(tfm, origin2DicomTransformFile)
tfm = self.applySlicer2AutoscoperTransform(tfm, segmentVolumeSize, segmentVolume.GetOrigin())
filename = f"{segmentName}.tra"
filename = os.path.join(outputDir, trackingSubDir, filename)
IO.writeTRA(filename, [tfm])

# update progress bar
progressCallback((idx + 1) / numSegments * 100)
Expand Down Expand Up @@ -1427,3 +1426,52 @@ def IsVolumeCentered(node: Union[slicer.vtkMRMLVolumeNode, slicer.vtkMRMLSequenc
if AutoscoperMLogic.IsSequenceVolume(node):
return AutoscoperMLogic.getItemInSequence(node, 0)[0].IsCentered()
return node.IsCentered()

@staticmethod
def applyOrigin2DicomTransform(transform: vtk.vtkMatrix4x4, Origin2DicomTransformFile: str) -> vtk.vtkMatrix4x4:
transformNode = slicer.vtkMRMLLinearTransformNode()
transformNode.SetMatrixTransformToParent(transform)
slicer.mrmlScene.AddNode(transformNode)
origin2DicomTransformNode = slicer.util.loadNodeFromFile(Origin2DicomTransformFile)
origin2DicomTransformNode.Inverse()
transformNode.SetAndObserveTransformNodeID(origin2DicomTransformNode.GetID())
transformNode.HardenTransform()
slicer.mrmlScene.RemoveNode(origin2DicomTransformNode)

transformMatrix = vtk.vtkMatrix4x4()
transformNode.GetMatrixTransformToParent(transformMatrix)
return transformMatrix

@staticmethod
def applySlicer2AutoscoperTransform(
transform: vtk.vtkMatrix4x4, volSize: list[float], origin: list[float]
) -> vtk.vtkMatrix4x4:
"""Utility function for converting a transform between the Slicer and Autoscoper coordinate systems."""
# https://github.com/BrownBiomechanics/Autoscoper/issues/280
transform.SetElement(1, 1, -transform.GetElement(1, 1)) # Flip Y
transform.SetElement(2, 2, -transform.GetElement(2, 2)) # Flip Z

# Add the volume origin
for i in range(3):
transform.SetElement(i, 3, transform.GetElement(i, 3) + origin[i])

transform.SetElement(0, 3, transform.GetElement(0, 3) - volSize[0]) # Offset X
return transform

@staticmethod
def vtkToNumpy(matrix: vtk.vtkMatrix4x4) -> NDArray:
"""Utility function for converting a 4x4 vtk matrix to a 4x4 numpy array."""
import numpy as np

array = np.empty((4, 4))
for i, j in product(range(4), range(4)):
array[i, j] = matrix.GetElement(i, j)
return array

@staticmethod
def numpyToVtk(array: NDArray) -> vtk.vtkMatrix4x4:
"""Utility function for converting a 4x4 numpy array to a 4x4 vtk matrix."""
matrix = vtk.vtkMatrix4x4()
for i, j in product(range(4), range(4)):
matrix.SetElement(i, j, array[i, j])
return matrix
63 changes: 11 additions & 52 deletions AutoscoperM/AutoscoperMLib/IO.py
Original file line number Diff line number Diff line change
Expand Up @@ -255,58 +255,6 @@ def writeTFMFile(filename: str, spacing: list[float], origin: list[float]):
slicer.mrmlScene.RemoveNode(transformNode)


def createTRAFile(
volName: str,
trialName: str,
outputDir: str,
trackingsubDir: str,
volSize: list[float],
Origin2DicomTransformFile: str,
transform: vtk.vtkMatrix4x4,
):
transformNode = slicer.vtkMRMLLinearTransformNode()
transformNode.SetMatrixTransformToParent(transform)
slicer.mrmlScene.AddNode(transformNode)

if Origin2DicomTransformFile is not None:
origin2DicomTransformNode = slicer.util.loadNodeFromFile(Origin2DicomTransformFile)
origin2DicomTransformNode.Inverse()
transformNode.SetAndObserveTransformNodeID(origin2DicomTransformNode.GetID())
transformNode.HardenTransform()
slicer.mrmlScene.RemoveNode(origin2DicomTransformNode)

filename = f"{trialName}_{volName}.tra" if trialName is not None else f"{volName}.tra"
filename = os.path.join(outputDir, trackingsubDir, filename)

if not os.path.exists(os.path.join(outputDir, trackingsubDir)):
os.mkdir(os.path.join(outputDir, trackingsubDir))

tfmMat = vtk.vtkMatrix4x4()
transformNode.GetMatrixTransformToParent(tfmMat)

writeTRA(filename, volSize, tfmMat)

slicer.mrmlScene.RemoveNode(transformNode)


def writeTRA(filename: str, volSize: list[float], transform: vtk.vtkMatrix4x4):
# Slicer 2 Autoscoper Transform
# https://github.com/BrownBiomechanics/Autoscoper/issues/280
transform.SetElement(1, 1, -transform.GetElement(1, 1)) # Flip Y
transform.SetElement(2, 2, -transform.GetElement(2, 2)) # Flip Z

transform.SetElement(0, 3, transform.GetElement(0, 3) - volSize[0]) # Offset X

# Write TRA
rowwise = []
for i in range(4): # Row
for j in range(4): # Col
rowwise.append(str(transform.GetElement(i, j)))

with open(filename, "w+") as f:
f.write(",".join(rowwise))


def writeTemporyFile(filename: str, data: vtk.vtkImageData) -> str:
"""
Writes a temporary file to the slicer temp directory
Expand Down Expand Up @@ -336,3 +284,14 @@ def removeTemporyFile(filename: str):

slicerTempDirectory = slicer.app.temporaryPath
os.remove(os.path.join(slicerTempDirectory, filename))


def writeTRA(fileName: str, transforms: list[vtk.vtkMatrix4x4]) -> None:
from itertools import product

rowWiseStrings = []
for transform in transforms:
rowWiseStrings.append([str(transform.GetElement(i, j)) for i, j in product(range(4), range(4))])
with open(fileName, "w") as traFile:
for row in rowWiseStrings:
traFile.write(",".join(row) + "\n")
3 changes: 2 additions & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ set(EXTENSION_CONTRIBUTORS "Anthony Lombardi (Kitware), Amy Morton (Brown Univer
set(EXTENSION_DESCRIPTION "SlicerAutoscoperM is an extension for 3D Slicer for 2D-3D image registration integrating with Autoscoper for image-based 3D motion tracking of skeletal structures.")
set(EXTENSION_ICONURL "https://raw.githubusercontent.com/BrownBiomechanics/SlicerAutoscoperM/main/SlicerAutoscoperM.png")
set(EXTENSION_SCREENSHOTURLS "https://github.com/BrownBiomechanics/SlicerAutoscoperM/releases/download/docs-resources/AutoscoperMainWindow.png")
set(EXTENSION_DEPENDS Sandbox SegmentEditorExtraEffects)
set(EXTENSION_DEPENDS Sandbox SegmentEditorExtraEffects SlicerElastix)
set(EXTENSION_BUILD_SUBDIRECTORY inner-build)

set(SUPERBUILD_TOPLEVEL_PROJECT inner)
Expand Down Expand Up @@ -44,6 +44,7 @@ add_subdirectory(AutoscoperM)
add_subdirectory(TrackingEvaluation)
add_subdirectory(CalculateDataIntensityDensity)
add_subdirectory(VirtualRadiographGeneration)
add_subdirectory(Tracking3D)
## NEXT_MODULE

#-----------------------------------------------------------------------------
Expand Down
34 changes: 34 additions & 0 deletions Tracking3D/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
#-----------------------------------------------------------------------------
set(MODULE_NAME Tracking3D)

#-----------------------------------------------------------------------------
set(MODULE_PYTHON_SCRIPTS
${MODULE_NAME}.py
${MODULE_NAME}Lib/__init__.py
${MODULE_NAME}Lib/TreeNode.py
)

set(MODULE_PYTHON_RESOURCES
Resources/Icons/${MODULE_NAME}.png
Resources/UI/${MODULE_NAME}.ui
Resources/ParameterFiles/rigid.txt
)

#-----------------------------------------------------------------------------
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()
Binary file added Tracking3D/Resources/Icons/Tracking3D.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
58 changes: 58 additions & 0 deletions Tracking3D/Resources/ParameterFiles/rigid.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
(FixedInternalImagePixelType "float")
(MovingInternalImagePixelType "float")

// **************** Main Components **************************

(Registration "MultiResolutionRegistration")
(Interpolator "BSplineInterpolator")
(ResampleInterpolator "FinalBSplineInterpolator")
(Resampler "DefaultResampler")
(FixedImagePyramid "FixedSmoothingImagePyramid")
(MovingImagePyramid "MovingSmoothingImagePyramid")

(Optimizer "AdaptiveStochasticGradientDescent")
(Transform "EulerTransform")
(Metric "AdvancedMattesMutualInformation")

// ***************** Transformation **************************

(AutomaticScalesEstimation "true")
(AutomaticTransformInitialization "true")
(HowToCombineTransforms "Compose")

//Save composite ITK transform
(ITKTransformOutputFileNameExtension "h5")
(WriteITKCompositeTransform "true")

// ******************* Similarity measure *********************

(NumberOfHistogramBins 64)
(ErodeMask "false")

// ******************** Multiresolution **********************

(NumberOfResolutions 3)
(ImagePyramidSchedule 8 8 2 4 4 1 1 1 1 )

// ******************* Optimizer ****************************

(MaximumNumberOfIterations 1000)
(MaximumStepLength 6.0)
(RequiredRatioOfValidSamples 0.01)

// **************** Image sampling **********************

(NumberOfSpatialSamples 2000)
(NewSamplesEveryIteration "true")
(ImageSampler "Random")

// ************* Interpolation and Resampling ****************

(BSplineInterpolationOrder 1)
(FinalBSplineInterpolationOrder 3)

(DefaultPixelValue 0)
(WriteResultImage "true")
(ResultImagePixelType "short")
(ResultImageFormat "mhd")

Loading

0 comments on commit b3a8fc1

Please sign in to comment.