search

SuperElastix / SimpleElastix

Multi-lingual medical image registration library
http://simpleelastix.github.io
Apache License 2.0
509 stars 149 forks source link

can't get deformation vector field #259

Open tiaplautz opened 5 years ago

tiaplautz commented 5 years ago

Hello,

When I run the suggested commands:

transformixImageFilter = sitk.TransformixImageFilter() transformixImageFilter.ComputeDeformationFieldOn() transformixImageFilter.SetTransformParameterMap(elastixImageFilter.GetTransformParameterMap()) transformixImageFilter.Execute() deformationField = transformixImageFilter.GetDeformationField()

it throws the runtime error:

RuntimeError Traceback (most recent call last)

in () 2 transformixImageFilter.ComputeDeformationFieldOn() 3 transformixImageFilter.SetTransformParameterMap(elastixImageFilter.GetTransformParameterMap()) ----> 4 transformixImageFilter.Execute() 5 deformationField = transformixImageFilter.GetDeformationField() C:\Users\tplautz\anaconda3\lib\site-packages\simpleitk-1.0.1rc1.dev345+g9dfa8-py3.7-win-amd64.egg\SimpleITK\SimpleITK.py in Execute(self) 10451 def Execute(self): 10452 """Execute(TransformixImageFilter self) -> Image""" > 10453 return _SimpleITK.TransformixImageFilter_Execute(self) 10454 10455 RuntimeError: Exception thrown in SimpleITK TransformixImageFilter_Execute: z:\simpleelastix\code\elastix\src\sitktransformiximagefilterimpl.cxx:109: sitk::ERROR: itk::ExceptionObject (000000506B5E81E0) Location: "unknown" File: z:\el\elastix\core\main\elxtransformixfilter.hxx Line: 218 Description: itk::ERROR: Self(000001B34E182890): Internal transformix error: See transformix log (use LogToConsoleOn() or LogToFileOn()) When I turn the LogToConsoleOn the output reads: ELASTIX version: 4.801 Command line options from ElastixBase: -out ./ -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 Calling all ReadFromFile()'s ... --------------- Exception --------------- itk::ExceptionObject (000000506B5E7890) Location: "unknown" File: z:\el\itk-prefix\include\itk-4.12\itkimagebase.hxx Line: 194 Description: itk::ERROR: Image(000001B34E052610): Bad direction, determinant is 0. Direction is 1 0 0 0 ----------------------------------------- Any advice on how to fix? Kindest Regards, Tia
kaspermarstal commented 5 years ago

Try to also set a moving image with the correct number of dimensions before executing. Let me know how it goes.

nwschurink commented 5 years ago

Hi, I am also having problems with the deformation field. Whenever I try to access the numpy data within the ITK vector image python crashes. 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. 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 parameter 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 pixel data works fine (e.g. deform.GetNumberOfComponentsPerPixel(), deform.GetSize() etc). However, when I use for example 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?

Also another peculiarity is that although I've set the output directory for transformix, no files are being written there. After running transformix through SimpleElastix the directory is still empty. When running the transformix command from command line itself everything goes fine and I get the resulting deformation field.

tiaplautz commented 5 years ago

@kaspermarstal thank you for your reply. I have added the line:

transformixImageFilter.SetMovingImage(movingImg)

Now the output reads:

Reading input image ... Reading input image took 0.000083 s Calling all ReadFromFile()'s ... Calling all ReadFromFile()'s took 0.016342 s Transforming points ... The transform is evaluated on all points. The result is a deformation field. Transforming points done, it took 0.22s 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 0.30s

but python throws the error: AttributeError: 'TransformixImageFilter' object has no attribute 'GetDeformationField'

-Tia

kaspermarstal commented 5 years ago

@tiaplautz It looks like you have an older version of SimpleElastix. Try to compile from scratch and install again.

kaspermarstal commented 5 years ago

@nwschurink Yes, it looks like the error is related to internal memory management. Will take a look. If you would be so kind as to open a separate issue?

tiaplautz commented 5 years ago

will do! Thanks very much!

MarHHM commented 5 years ago

@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 :)

tychovdo commented 4 years ago

I also struggle to access the pixel data of the deformation field after calling transformixImageFilter.GetDeformationField(), resulting in a crash without raising an error. The issue seems to be similar to the one of @nwschurink.

@MarHHM @kaspermarstal @nwschurink Any ideas how to fix this issue?

Thanks in advance.

AmericaBG commented 4 years ago

Hi!! And finally...Had you find how to fix the error?😅

Thank you very much!

nwschurink commented 4 years ago

Hi!! And finally...Had you find how to fix the error?😅

Thank you very much!

Hi @SouthAmericaB, I fixed it by building superelastix from the dev branch (the bug is fixed as mentioned here https://github.com/SuperElastix/SimpleElastix/issues/261#issuecomment-593807204 and https://github.com/SuperElastix/SimpleElastix/issues/261#issuecomment-460596708) Another (crude) sollution is to use subprocess to call elastix and transformix from command line and load the result using ITK.

AmericaBG commented 4 years ago

Ey! Thank you very much for your reply! I built superelastix from the dev branch a days ago. I'm using Anaconda and Spyder environment for Python code.

My code: `from future import print_function import SimpleITK as sitk import sys import os

variable='G:/Datos_STIR_SagT2/STIR/1.nii' fija='G:/Datos_STIR_SagT2/SagitalT2/1_Sag.nii'

elastixImageFilter = sitk.ElastixImageFilter() elastixImageFilter.SetFixedImage(sitk.ReadImage(fija)) elastixImageFilter.SetMovingImage(sitk.ReadImage(variable)) elastixImageFilter.SetParameterMap(sitk.GetDefaultParameterMap("affine")) elastixImageFilter.Execute()

transformixImageFilter = sitk.TransformixImageFilter() transformixImageFilter.SetTransformParameterMap(elastixImageFilter.GetTransformParameterMap()) transformixImageFilter.ComputeDeformationFieldOn() transformixImageFilter.Execute() deformationField = transformixImageFilter.GetDeformationField()

parameterMapVector = elastixImageFilter.GetParameterMap()

sitk.WriteImage(elastixImageFilter.GetResultImage(),'G:/Datos_STIR_SagT2/STIR/fds.nii')

sitk.WriteImage(deformationField,'G:/Datos_STIR_SagT2/STIR/de.nii') ¿????????????

` My problem is: the code run correctly but I can't see (or better, I don't know how) the "deformationField variable" .

Could you explain me how do you do to see this result???

Thank you very very much!

nwschurink commented 4 years ago

@SouthAmericaB , I think you mean you want to visualize the deformationfield?

You can use software such as 3D slicer. Import the deformationfield as a 'Transform' image. After importing it you can visualize the transform as described here.

AmericaBG commented 4 years ago

Sorry, I think I haven't explained very well 🤦‍♀ I want to extract information of the deformationfield, either as an image, as a matrix, the pixel values... But I don't find where this variable is saved after my code have run. I have searched in the working directory, in the SimpleElastix directory...but there's nothing!

nwschurink commented 4 years ago

You can save the resulting deformation field as a file using the code you commented out:

sitk.WriteImage(deformationField,'G:/Datos_STIR_SagT2/STIR/de.nii')

Or another option is that you use the command transformixImageFilter.SetOutputDirectory("./path/to/output/directory") , before you run transformixImageFilter.Execute(). It should then save all files to the path that you give. After that you can process the files you would normally do.

Or , if you want to get the pixel values as a numpy matrix in python directly you can probably do something like :

deformationField_array = sitk.GetArrayFromImage(deformationField)

deformationField_array should then contain your numpy array.

AmericaBG commented 4 years ago

Thank you ver much! I didn't know some of that information (for example, to get values as a numpy array)!

My error was that I wasn't assigning a value to transformixImageFilter.SetMovingImage() that's why I got no value after transformixImageFilter.Execute() (no deformationField, no transformixImageFilter.GetResultImage()...)

Now, the code performs correctly.

Thank you very much for your time and attention :)

ctorti commented 4 years ago

@nwschurink and @SouthAmericaB ,

Thanks to Niels for your suggestions to convert the deformation field to a numpy array, or to write the image to disk. I tried both and unfortunately in both cases it killed my kernel.

Any ideas why?

My system details:

Anaconda Jupyter 6.0.1, IPython 7.8.0, Python 3.7.4, SimpleITK 1.2.0rc2.dev1162 (ITK 4.13), Win10 64-bit

And @SouthAmericaB: I'm interesed to see that by setting the moving image you got the expected output after running Transformix. Unfortunately for me that hasn't made any difference. May I ask what IDE are you using?

Has anyone managed to get the output to work when using the same/similar system configuration as I have shown above?

Aside: I'm not sure that this is relevant, but despite setting .LogToConsoleOn() I don't get any output to my console. And as far as and .LogToFileOn(), the log file elastix.log is only written/updated once per kernel session for some reason.

I've seen some other threads discussing these issues in Anaconda but have yet to find a workable solution. E.g. 

import sys
sys.stdout.flush()

didn't work for me.

So I'm stuck in a rather difficult situation whereby:

Does anyone have any ideas?

Here's my code if interested - first the registration section:

# 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']
# Print the parameters:
for keys,values in ParamMap.items():
    print(keys, '=', values)
# Set the parameter map:
ElastixImFilt.SetParameterMap(ParamMap)

# Register the 3D images:
ElastixImFilt.Execute()
# Get the registered image:
RegIm = ElastixImFilt.GetResultImage()

I'm happy with the above - the registration looks fine.

And now, my attempts 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()
#TransformixIm = TransformixImFilt.Execute()

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

DeformationField

Output:

<SimpleITK.SimpleITK.Image; proxy of <Swig Object of type 'std::vector< itk::simple::Image >::value_type *' at 0x0000026793233750> >

Thanks!

ApoorvaDA commented 3 years ago

@SouthAmericaB , I think you mean you want to visualize the deformationfield?

You can use software such as 3D slicer. Import the deformationfield as a 'Transform' image. After importing it you can visualize the transform as described here.

In this case, how can we save the deformation field as a 'tfm' file? Can 3D Slicer read a nifty file as transform file? Thanks, Apoorva.D.A

lassoan commented 3 years ago

Nifti fine format had many issues but if you save the field as a nrrd file then 3D Slicer can read it. You can also run Elastix in 3D Slicer with a convenient GUI or by a few lines of Python script by installing SlicerElastix extension. It uses a pre-built Elastix that is compatible with 3D Slicer's Python environment, so you won't have crashes and you don't need to build Elastix yourself.