CoBrALab / optimized_antsMultivariateTemplateConstruction

A re-implementation of antsMultivariateTemplateConstruction2.sh using optimized image pyramid scale-space and qbatch support
Other
22 stars 8 forks source link

Introduce scaling for affine transforms during modelbuild #13

Closed gdevenyi closed 1 year ago

gdevenyi commented 3 years ago

Implementation in python, depends on sitk and vtk

#!/usr/bin/env python

# https://vtk.org/doc/nightly/html/classvtkTransformInterpolator.html
# https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1Transform.html
# https://vtk.org/doc/nightly/html/classvtkMatrix4x4.html
# https://www.paraview.org/Bug/view.php?id=10102#c37133

import vtk
import SimpleITK as sitk

# Input transform
input_itk_transform = sitk.ReadTransform("trans0_GenericAffine.mat")
output_itk_transform = sitk.AffineTransform(3)
# Dump the matrix 4x3 matrix
input_itk_transform_paramters = input_itk_transform.GetParameters()

# Create VTK input/output transform, and an identity for the interpolation
input_vtk_transform = vtk.vtkTransform()
identity_vtk_transform = vtk.vtkTransform()
output_vtk_transform = vtk.vtkTransform()

# Setup the VTK transform by reconstructing from the ITK matrix parameters
input_vtk_transform.SetMatrix(
    input_itk_transform_paramters[0:3]
    + input_itk_transform_paramters[9:10]
    + input_itk_transform_paramters[3:6]
    + input_itk_transform_paramters[10:11]
    + input_itk_transform_paramters[6:9]
    + input_itk_transform_paramters[11:12]
    + (0, 0, 0, 1)
)

# Create an interpolator
vtk_transform_interpolator = vtk.vtkTransformInterpolator()

# Build an interpolation stack, identity transform and the input transform
vtk_transform_interpolator.AddTransform(0, identity_vtk_transform)
vtk_transform_interpolator.AddTransform(1, input_vtk_transform)

for scale in [0.1,0.25,0.5,0.75,1]:
    # Generate a transform a fractional step between identity and the input transform
    vtk_transform_interpolator.InterpolateTransform(scale, output_vtk_transform)

    output_itk_transform.SetParameters(
        (
            output_vtk_transform.GetMatrix().GetElement(0, 0),
            output_vtk_transform.GetMatrix().GetElement(0, 1),
            output_vtk_transform.GetMatrix().GetElement(0, 2),
            output_vtk_transform.GetMatrix().GetElement(1, 0),
            output_vtk_transform.GetMatrix().GetElement(1, 1),
            output_vtk_transform.GetMatrix().GetElement(1, 2),
            output_vtk_transform.GetMatrix().GetElement(2, 0),
            output_vtk_transform.GetMatrix().GetElement(2, 1),
            output_vtk_transform.GetMatrix().GetElement(2, 2),
            output_vtk_transform.GetMatrix().GetElement(0, 3),
            output_vtk_transform.GetMatrix().GetElement(1, 3),
            output_vtk_transform.GetMatrix().GetElement(2, 3)
        )
    )

    sitk.WriteTransform(output_itk_transform, f"output{str(scale)}.mat")