Can't access or save data array in transformix deformation vector field

[nwschurink]

Hi, I am posting this a new issue as requested by @kaspermarstal . I am having problem's with the transformix module of SimpleElastix. More specifically, whenever I try to access the numpy data within the ITK deformation vector field image python crashes without raising an error. Another peculiarity is that Transformix is not saving any files to the defined output folder, wheras Elastix is saving files to the correct folder.

I use the following code:

selx = sitk.ElastixImageFilter()
selx.SetParameterMap(selx.ReadParameterFile("./Parameters_BSpline.txt"))
selx.SetOutputDirectory("./output")
fixedImage = sitk.ReadImage("./T2_TRA_image.nii.gz")
fixedMask = sitk.ReadImage("./TUMOR_ROI.nii.gz")
# There is a small offset (1e-10) in the direction cosines which lets elastix fail. Therefore we resample
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixedImage)
fixedMask = resampler.Execute(fixedMask)
selx.SetFixedImage(fixedImage)
selx.SetFixedMask(sitk.Cast(fixedMask, sitk.sitkUInt8))
movingImage = sitk.ReadImage("./CT_to_T2.nii.gz")
selx.SetMovingImage(movingImage)
selx.Execute()
transformedImage = selx.GetResultImage()

strx = sitk.TransformixImageFilter()
strx.SetMovingImage(movingImage)
strx.ComputeDeformationFieldOn()
strx.LogToConsoleOn()
strx.SetOutputDirectory("./deform_output")
strx.Execute()

deform = strx.GetDeformationField()

Up until here everything seems to work fine. The output of transformix to the console is:

ELASTIX version: 4.900
Command line options from ElastixBase:
-out      ./deform_output/
-threads  unspecified, so all available threads are used
-def      all
-jac      unspecified, so no det(dT/dx) computed
-jacmat   unspecified, so no dT/dx computed
Reading input image ...
  Reading input image took 0.000022 s
Calling all ReadFromFile()'s ...
WARNING: The parameter "UseBinaryFormatForTransformationParameters", requested at entry number 0, does not exist at all.
  The default value "false" is used instead.
  Calling all ReadFromFile()'s took 0.016059 s
Transforming points ...
  The transform is evaluated on all points. The result is a deformation field.
  Transforming points done, it took 1.43s
Compute determinant of spatial Jacobian ...
  The command-line option "-jac" is not used, so no det(dT/dx) computed.
  Computing determinant of spatial Jacobian done, it took 0.00s
Compute spatial Jacobian (full matrix) ...
  The command-line option "-jacmat" is not used, so no dT/dx computed.
  Computing spatial Jacobian done, it took 0.00s
Resampling image and writing to disk ...
  Resampling took 1.65s
Out[6]: <SimpleITK.SimpleITK.Image; proxy of <Swig Object of type 'std::vector< itk::simple::Image >::value_type *' at 0x0000000005C3C8A0> >

Ok so just to get some information about the deform variable that I defined, I printed it's info which is:

VectorImage (0000000006620DF0)
  RTTI typeinfo:   class itk::VectorImage<float,3>
  Reference Count: 2
  Modified Time: 7180
  Debug: Off
  Object Name: 
  Observers: 
    none
  Source: (none)
  Source output name: (none)
  Release Data: Off
  Data Released: False
  Global Release Data: Off
  PipelineMTime: 0
  UpdateMTime: 0
  RealTimeStamp: 0 seconds 
  LargestPossibleRegion: 
    Dimension: 3
    Index: [0, 0, 0]
    Size: [256, 256, 35]
  BufferedRegion: 
    Dimension: 3
    Index: [0, 0, 0]
    Size: [256, 256, 35]
  RequestedRegion: 
    Dimension: 3
    Index: [0, 0, 0]
    Size: [256, 256, 35]
  Spacing: [0.78125, 0.78125, 5]
  Origin: [-103.11, 142.614, -87.0026]
  Direction: 
0.999869 -0.0017434 0.0160782
-0.00246245 -0.998992 0.0448114
-0.0159839 0.0448451 0.998866
  IndexToPointMatrix: 
0.781148 -0.00136203 0.080391
-0.00192379 -0.780462 0.224057
-0.0124874 0.0350352 4.99433
  PointToIndexMatrix: 
1.27983 -0.00315194 -0.0204594
-0.00223155 -1.27871 0.0574018
0.00321565 0.00896228 0.199773
  Inverse Direction: 
0.999869 -0.00246245 -0.0159839
-0.0017434 -0.998993 0.0448452
0.0160782 0.0448114 0.998866
  VectorLength: 3
  PixelContainer: 
    ImportImageContainer (0000000004190530)
      RTTI typeinfo:   class itk::ImportImageContainer<unsigned __int64,float>
      Reference Count: 1
      Modified Time: 7173
      Debug: Off
      Object Name: 
      Observers: 
        none
      Pointer: 000000000B9E0040
      Container manages memory: false
      Size: 6881280
      Capacity: 6881280

Using commands that don't try to access the numpy pixel data works fine (e.g. deform.GetNumberOfComponentsPerPixel(), deform.GetSize() etc). However, when I use deform.GetPixel((1,1,1)) or sitk.GetArrayFromImage(deform), or when I try to write the deformation field to a file using sitk.WriteImage(deform, 'deform.mhd') python crashes without raising an error.

The exact version of SimpleElastix that I have installed is simpleitk-1.1.0.dev362+gb3783-py3.6-win-amd64 according to conda. Could it be that this is some bug in the C code?

When running the transformix command from command line itself everything goes fine and I get the resulting deformation field.

[kaspermarstal]

It should now work on the develop branch. New tests pass on CI. Give it a spin and let me know how it goes.

[MarHHM]

@kaspermarstal I got exactly the same issue as @nwschurink .. I can't write the deformation field I get through TransformixImageFilter.GetDeformationField() with sitk.WriteImage().The generated .raw is basically empty, although the generated .mhd seems to have the correct metadata.

Also, I can't access the pixel data within the deformation field as well..

I'd be happy if there are any updates regarding this issue :)

Thanks for the great effort :)

[nwschurink]

@MarHHM It is a while ago that I have encountered it, but I think it was solved on the dev branch as is mentioned by @kaspermarstal.
Building SuperElastix from dev should probably fix your problem.

Another (crude) solution (which was a workaround I resorted to at first) would be to run the transformix command from command line (using eg subprocess) and load the deformation field using SimpleITK. It's a messy workaround though, so probably better to just use the dev build.

[MarHHM]

@MarHHM It is a while ago that I have encountered it, but I think it was solved on the dev branch as is mentioned by @kaspermarstal.
Building SuperElastix from dev should probably fix your problem.

Another (crude) solution (which was a workaround I resorted to at first) would be to run the transformix command from command line (using eg subprocess) and load the deformation field using SimpleITK. It's a messy workaround though, so probably better to just use the dev build.

Thanks for ur feedback, @nwschurink

A while ago, I resorted to using the command line versions (v5.0) of both elastix & transformix (calling them in python through subprocess.run() ) that I build separately without the Superbuild pipeline of sitk (as the version that sitk uses till now is still 4.9 which has performance issues & also ElastixImageFilter doesn't provide a way to supply some command line args (e.g. "-labels" used for (Transform "MultiBSplineTransformWithNormal") that i was interested in)).

I use sitk currently only for pre & postprocessing (i.e. ITK functionality, but not elastix) as it's still super useful..

Cheers,,
M

[ctorti]

Hi all,

I've gone through quite a painful process to finally get a deformation file exported. But unfortunately the entries are all zero.

I initially performed a registration in Anaconda Jupyter:

# Initiate ElastixImageFilter:
ElastixImFilt = sitk.ElastixImageFilter()
ElastixImFilt.LogToConsoleOn() # <-- no output
ElastixImFilt.LogToFileOn() # <-- output's elastix.log ONLY ONCE per kernel

# Define the fixed and moving images:
ElastixImFilt.SetFixedImage(FixedIm)
ElastixImFilt.SetMovingImage(MovingIm)

# Get the default parameter map template:
ParamMap = sitk.GetDefaultParameterMap('affine')

# Set registration parameters:
ParamMap['AutomaticTransformInitialization'] = ['true']
ParamMap['AutomaticTransformInitializationMethod'] = ['GeometricalCenter']
ParamMap['WriteIterationInfo'] = ['true']
ParamMap['MaximumNumberOfIterations'] = ['512']
ParamMap['UseDirectionCosines'] = ['true']

# Set the parameter map:
ElastixImFilt.SetParameterMap(ParamMap)

# Register the 3D images:
ElastixImFilt.Execute()

# Get the registered image:
RegIm = ElastixImFilt.GetResultImage()

and the registration looked fine. So I moved onto trying to get the deformation field:

# Get the current working directory:
CurrentWorkingDir = os.getcwd()

# Get the transform parameter map:
#TransformParameterMap = ElastixImFilt.GetParameterMap()

# Initiate TransformixImageFilter:
TransformixImFilt = sitk.TransformixImageFilter()
TransformixImFilt.LogToConsoleOn() # <-- no output
TransformixImFilt.LogToFileOn() 
TransformixImFilt.SetOutputDirectory(CurrentWorkingDir) # <-- makes no difference, still no output

# Set up for computing deformation field:
TransformixImFilt.ComputeDeformationFieldOn() # <-- no output

# Set up other computations to see if they produce output:
TransformixImFilt.ComputeSpatialJacobianOn() # <-- produces 'spatialJacobian.nii'
TransformixImFilt.ComputeDeterminantOfSpatialJacobianOn() # <-- produces 'fullSpatialJacobian.nii'

# Set the transform parameter map:
TransformixImFilt.SetTransformParameterMap(ElastixImFilt.GetTransformParameterMap())

# Need to explicitely tell elastix that the image is 3D since by default it will
# only transform to 2D:
TransformixImFilt.SetMovingImage(MovingIm)

# Transform the contour points:
TransformixImFilt.Execute()

# Get the deformation field:
DeformationField = TransformixImFilt.GetDeformationField()

but the deformation field was not exported (the Jacobian and its determinant was exported to .nii files).

Issue #259 looked promising so I tried:

DeformationFieldNda = sitk.GetArrayFromImage(DeformationField)

but it killed the kernel without producing any errors. So did:

sitk.WriteImage(DeformationField, 'defField.nii')

Then I copied the above code and saved it into a .py file and ran it in Anaconda Spyder. Spyder produced output (warnings, parameters, metrics results, etc) in the console that Jupyter did not, but still not deformationField.nii file.

Then I repeated the above in PyCharm and it worked! Or so I thought it did...

Running the script in PyCharm created the defField.nii file only once! Although the Jacobian .nii files and other .log and .txt files were updated, the deformation field file did not. I tried again after restarting PyCharm but there was no change.

I moved on. I read in the defField.nii file I generated once only and converted to a numpy array:

DefField = sitk.ReadImage(DefFieldFpath)
DefFieldNda = sitk.GetArrayFromImage(DefField)

but found that all entries are 0. It's quite a defeating result. Does anyone have any idea what's going on?

As far as I'm aware I'm running the latest build of SimpleElastix (1.2.0rc2.dev1162-g2a79d).

Thanks.