Extracting displacement fields from rigid registration

Hello!

Continuing the discussion from 3D image registration with Python and without SimpleITK:

Once I have a registration, what is the best way to extract the deformation field (a.k.a. displacement field) with ITK Python? The idea is to avoid reprocessing and have shorter code.

  1. is converting my image to points still necessary (see this).
    • is this also the case once the scaling, origin and direction are set for the image (as shown in Sec. 4.14 Defining Origin and Spacing of the Software Guide)?
    • does one need to iterate over the whole image (to use points)?
  2. is a composite transform (permalink) a better alternative? (As I am looking for a Python ITK solution, I would translate the SimpleITK syntax).
  3. is there a way to get the displacement (deformation) field from the registration or transform directly? (as in: a method of the Transform class) [Edit]. itk.TransformToDisplacementFieldFilter tells me that AttributeError: module 'itk' has no attribute 'TransformToDisplacementFieldFilter'
  4. is there a full example which shows how to extract the displacements from a rigid registration?

Thank you, very much!

TransformToDisplacementFieldFilter is wrapped for Python. If you compiled ITK from source, make sure that ITK_WRAP_vector_float is enabled and ITK_WRAP_VECTOR_COMPONENTS contains 3. A Python example which touches a deformation field might be useful from syntax perspective.

1 Like

Hi there, Dženan. Thank you. My CMakeCache.txt shows:

//Number of vector components available separated by semicolons
// (;)
ITK_WRAP_VECTOR_COMPONENTS:STRING=2;3;4
...
//Wrap vector float type
ITK_WRAP_vector_float:BOOL=ON

and my compilation instructions include

-DITK_WRAP_IMAGE_DIMS:STRING="2;3;4"

Do you think that there would be a difference if I recompile like this?:

-DITK_WRAP_VECTOR_COMPONENTS:STRING="2;3;4" -DITK_WRAP_vector_float:BOOL=ON

I don’t think there would be a difference. My cache has the same content, I think those are default options. There is no problem for me:

>>> import itk
>>> itk.TransformToDisplacementFieldFilter
<itkTemplate itk::TransformToDisplacementFieldFilter>
>>>

Thank you again, Dženan. As I am looking into my makepkg.conf (makepkg is the building tool of Arch-based systems), I see that there is an option to disable and clear static libraries. I wonder if that has to do with it. I will disable the option, compile and report back. Thanks.

I will create a question to solve my import situation and come back to this once I am done

1 Like

Hi! I’m back :slight_smile: .

I’ve got a little step forward:

displ_filt = itk.TransformToDisplacementFieldFilter[itk.Image[itk.Vector[itk.F,3],3], itk.D].New()

In the example which is linked above, the deformation field is provided as an input (deformationField = fieldReader.GetOutput()). I already tried to look into the documentation mentioned by @matt.mccormick over here, and with help(displ_filt). Right now, I am unsure of how to get the deformation field from the registration or movingInitialTransform. Any tips? Thank you in advance.

It’s nice to be back.

I think you only need these 2 lines more:

displ_filt.Update()
deformationField = displ_filt.GetOutput()
1 Like

Thank you, Dženan. I initially got this error:

RuntimeError: /pkg/insight-toolkit/src/insight-toolkit/Modules/Core/Common/src/itkProcessObject.cxx:1339:
ITK ERROR: TransformToDisplacementFieldFilter(0x55578b3cac50): Input Transform is required but not set.

Then, I used this instead:

displ_filt.SetTransform(registr_transform)
displ_filt.Update()
displ_field_out = displ_filt.GetOutput()

To save to a file, I naïvely used this:

displ_writer = itk.ImageFileWriter.New(
    Input=displ_field_out, FileName=displ_field_vtk)
displ_writer.Update()

that produces:

RuntimeError: /pkg/insight-toolkit/src/insight-toolkit/Modules/IO/ImageBase/include/itkImageFileWriter.hxx:245:
ITK ERROR: ImageFileWriter(0x555788553d50): Largest possible region does not fully contain requested paste IO regionPaste IO region: ImageIORegion (0x7ffd35280600)
  Dimension: 3
  Index: 0 0 0
  Size: 0 0 0
Largest possible region: ImageRegion (0x7ffd35280760)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [0, 0, 0]

I tried other approaches, but I’m sure that expanding on my mistakes is of no use :stuck_out_tongue: . How do I save the results? (is what I did to get them the right way?). Thank you again :slight_smile: .

This approach might need explicit type specification, such as displ_writer = itk.ImageFileWriter[itk.Image[itk.Vector[itk.D,3],3]].New(.

Simpler, recommended way is itk.imwrite(displ_field_out, displ_field_vtk)

1 Like

Thanks again. (I changed itk.D to itk.F) I am still getting the same error (with different memory allocation, possibly).

>>> itk.imwrite(displ_field_out, displ_field_vtk)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.10/site-packages/itk/support/extras.py", line 910, in imwrite
    writer.Update()
RuntimeError: /home/edgar/Documentos/Linux/Parabola/pkg/insight-toolkit/src/insight-toolkit/Modules/IO/ImageBase/include/itkImageFileWriter.hxx:245:
ITK ERROR: ImageFileWriter(0x55b1c49b2570): Largest possible region does not fully contain requested paste IO regionPaste IO region: ImageIORegion (0x7fff1c528160)
  Dimension: 3
  Index: 0 0 0 
  Size: 0 0 0 
Largest possible region: ImageRegion (0x7fff1c5282c0)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [0, 0, 0]

>>> displ_writer = itk.ImageFileWriter[itk.Image[itk.Vector[itk.F, 3], 3]].New(Input=displ_field_out, FileName=displ_field_vtk)
>>> displ_writer.Update()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
RuntimeError: /home/edgar/Documentos/Linux/Parabola/pkg/insight-toolkit/src/insight-toolkit/Modules/IO/ImageBase/include/itkImageFileWriter.hxx:245:
ITK ERROR: ImageFileWriter(0x55b1c4b2ef90): Largest possible region does not fully contain requested paste IO regionPaste IO region: ImageIORegion (0x7fff1c528330)
  Dimension: 3
  Index: 0 0 0 
  Size: 0 0 0 
Largest possible region: ImageRegion (0x7fff1c528490)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [0, 0, 0]

Oh, you need to specify the size etc of the desired displacement field. Use TransformToDisplacementFieldFilter.SetReferenceImage() plus UseReferenceImageOn or SetOutputOrigin, SetOutputSpacing etc. See docs for SetReferenceImage.

1 Like

Great, thanks. Paraphrasing from the docs:

This method can be used to specify an image from which to copy the [output information (SetOutputSpacing, Origin, and Direction)]. UseReferenceImageOn must be set

So I did:

displ_filt = itk.TransformToDisplacementFieldFilter[
    itk.Image[itk.Vector[itk.F, 3], 3], itk.D].New()
displ_filt.SetTransform(registr_transform)
displ_filt.SetReferenceImage(fixed_image)
displ_filt.UseReferenceImageOn()
displ_filt.Update()
displ_field_out = displ_filt.GetOutput()

and

itk.imwrite(displ_field_out, displ_field_vtk)

It runs well now :slight_smile: .

The (VTK) output shows up in ParaView as a block with the size of my fixed_image, completely full of little glyphs in the same direction. This makes me think that I will have to do this with a moving window (an area smaller than the full image) whose position I will move iteratively to get the displacements at each point, correct? This is possibly what @matt.mccormick was trying to say:

1 Like

If you need it only in a subset of the image, it would take less time to produce it and less memory to store it.

Dear Dženan,

I would like to have the displacements (deformation) at each point.

With my moving smiles (the top-most is the reference; the one at the bottom is moved to the right; the difference between them is shown to the left)

, I was expecting to get zero deformation in the gray areas. With the results from TransformToDisplacementFieldFilter (by means of itk.imwrite(displ_field_out, displ_field_vtk)), I get a uniform field of vectors . I am assuming that these vectors are the values of deformation.

Since the field is uniform, what I am assuming is that the overall displacement of the smile is applied in the whole region. That is, having used fixed_image (which is the full original image) results in the overall value being applied everywhere.

Therefore, I would think that I need to calculate the deformation at each point… but isn’t that the goal of TransformToDisplacementFieldFilter? (to get the deformation at each point?). Thus, I must be doing something wrong.

1 Like

Displacement field provides translation vector at each point. But if you construct it by computing rigid transform for each of those points, you get a rather uniform DF. So if you have a rigid transform, it is best to keep it in some specialized format (a few numbers specifying translation, rotation and center of rotation).

What you might be looking for is deformable registration. But even with deformable registration, most algorithms would produce movement in the gray areas between the eyes. Why? Because if they didn’t, it would imply huge deformations around the eyes and mouth and no deformation for the rest of the face. Put another way: if your head is translated, each point of your forehead is translated too, not just the eyes.

2 Likes

As always, thank you very much. This will certainly help me.

1 Like