# Split deformation vector field into two deformation vector fields and apply to image in order

**URL:** https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148
**Category:** Algorithms
**Created:** [August 15, 2023, 6:09am UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148 "2023-08-15T06:09:59Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![Xinyue](https://discourse.itk.org/user_avatar/discourse.itk.org/xinyue/32/3339_2.png) [@Xinyue](https://discourse.itk.org/u/Xinyue)
#### Post date: [August 15, 2023, 6:09am UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/1 "2023-08-15T06:09:59Z")

</div>

Hi,

I am trying to split the dvf into two dvfs and apply in order.

However the result of deforming with one dvf is not the same as deforming with the split two dvfs.

I am wondering why is happen and how to fix it.

here is my code:

> def deform\_img(img,dvf):
> 
> ```
> resampler = sitk.ResampleImageFilter()
> resampler.SetReferenceImage(img)
> dis_tx = sitk.DisplacementFieldTransform(sitk.Cast(dvf, sitk.sitkVectorFloat64))
> resampler.SetTransform(dis_tx)
> new_rigid_post = resampler.Execute(img)
> return new_rigid_post
> 
> ```
> 
> def test10():
> 
> ```
> post = np.zeros((128, 256, 256))
> post[40:90, 40:90, 40:90] = 255
> itk_post = sitk.GetImageFromArray(post)
> 
> pre = np.zeros((128, 256, 256))
> pre[30:80, 20:70, 10:60] = 255
> itk_pre = sitk.GetImageFromArray(pre)
> 
> dvf = np.zeros((128, 256, 256, 3))
> dvf[:, :, :, 0] = 30.2
> dvf[:, :, :, 1] = 20.2
> dvf[:, :, :, 2] = 10.2
> itk_dvf = sitk.GetImageFromArray(dvf)
> new_post=deform_img(itk_post,itk_dvf)
> new_post_img = sitk.GetArrayFromImage(new_post)
> 
> dvf1 = np.zeros((128, 256, 256, 3))
> dvf1[:, :, :, 0] = 20.0
> dvf1[:, :, :, 1] = 15.0
> dvf1[:, :, :, 2] = 5.0
> itk_dvf1 = sitk.GetImageFromArray(dvf1)
> part1_post = deform_img(itk_post, itk_dvf1)
> 
> dvf2 = np.zeros((128, 256, 256, 3))
> dvf2[:, :, :, 0] = 10.2
> dvf2[:, :, :, 1] = 5.2
> dvf2[:, :, :, 2] = 5.2
> itk_dvf2 = sitk.GetImageFromArray(dvf2)
> part2_post = deform_img(part1_post, itk_dvf2)
> part2_post_img = sitk.GetArrayFromImage(part2_post)
> print(np.array_equal(new_post_img,part2_post_img))
> 
> ```

The print(np.array\_equal(new\_post\_img,part2\_post\_img)) returns false.

---

<div class="post-metadata">

### Author: ![blowekamp](https://discourse.itk.org/user_avatar/discourse.itk.org/blowekamp/32/79_2.png) [@blowekamp](https://discourse.itk.org/u/blowekamp)
#### Post date: [August 15, 2023, 1:20pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/2 "2023-08-15T13:20:44Z")

</div>

That is correct based on the mathematical properties of displacement fields. I don’t believe this is covered in the ITK software guide.

~~I believe dvf2 would need to be transformed/resampled by inverse of dvf1.~~ The math of compositing deformations fields needs further thought.

Hopefully some one else can recommend an appropriate text to review the Mathematica foundations of displacements fields. @matt.mccormick @dzenanz @gdevenyi @zivy

---

<div class="post-metadata">

### Author: ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [August 15, 2023, 1:51pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/3 "2023-08-15T13:51:03Z")

</div>

Hello @Xinyue,

The difference you are observing is not specific to displacement fields. You are experiencing the discreetness of computation vs. the continuous world. This is a feature of resampling multiple times as compared to compositing transformations and resampling once (still suffers from the issue but minimizes its effects). To see that the results are equal, equality is always up to a certain \epsilon, add these two lines to the end of `test10()`:

```auto
  # Show all the locations where the values differ and
  # print the sum of absolute differences
  sitk.Show(((part2_post - new_post)!=0)*255, "diff")
  print(np.abs(new_post_img - part2_post_img).sum())

```

You will see that the differences are at the boundary of the box object, due the resampling, and that these differences are very small (sum of absolute differences is 4.6e-09).

Finally, if you are using deep learning, everybody is, MONAI added what they refer to as [lazy resampling](https://docs.monai.io/en/stable/lazy_resampling.html) composing spatial transformations which reduces the resampling artifacts associated with consecutive resampling (which is what the provided code does).

---

<div class="post-metadata">

### Author: ![Xinyue](https://discourse.itk.org/user_avatar/discourse.itk.org/xinyue/32/3339_2.png) [@Xinyue](https://discourse.itk.org/u/Xinyue)
#### Post date: [August 15, 2023, 2:46pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/4 "2023-08-15T14:46:28Z")

</div>

Hi @zivy,

Thank you for the explanation.

I think MONAI would be helpful when people try to resample multiple times and get one output.

I am working on image registration and my data need to be “rigid registration” between the pre and post image. I have the deformation vector field (dvf in my sample code) to move the post to look like pre, but I want to split this deformation vector field into two vector fields (dvf1 and dvf2), the first dvf moves all the voxels on post to same distance and same direction (which is like “rigid registration”), and the second dvf will deform the new post to my target. dvf = first dvf +second dvf.

I work on my deep learning model, the old input would be the old post images and the pre images, and out put is the dvf, but the new input would be the new post and the pre images.

The hard part is I need the output image after applying the dvf1. I am wondering if MONAI could help this issue since it reduce the multiple resample affect by only execute once.

Thank you,  
Xinyue

---

<div class="post-metadata">

### Author: ![gdevenyi](https://discourse.itk.org/user_avatar/discourse.itk.org/gdevenyi/32/355_2.png) [@gdevenyi](https://discourse.itk.org/u/gdevenyi)
#### Post date: [August 15, 2023, 4:33pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/5 "2023-08-15T16:33:22Z")

</div>

> [@Xinyue](#):
>
> I am working on image registration and my data need to be “rigid registration”

Why not compute an actual rigid transformation first?

> I am wondering if MONAI could help this issue since it reduce the multiple resample affect by only execute once.

You can reduce resampling error by using better interpolation techniques, but fundamentally sampling into a new space will always add some computational error.

---

<div class="post-metadata">

### Author: ![Xinyue](https://discourse.itk.org/user_avatar/discourse.itk.org/xinyue/32/3339_2.png) [@Xinyue](https://discourse.itk.org/u/Xinyue)
#### Post date: [August 15, 2023, 9:00pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/6 "2023-08-15T21:00:46Z")

</div>

Hi @gdevenyi ,

Thank you for your reply.

I am trying to use a simple way to move the image slightly, so the value of my deformation vector field would not be too large, and the data distribution would look better.

I don’t know if I use a rigid registration step in my data preprocessing will slow down the whole pipeline cause I want to finish the data preprocessing by real time.

I was trying to use 3D Slicer to see if my assumption (dvf = dvf1+dvf2) is correct, and it seems the 3D slicer won’t have any issue with deforming two times, but I am not sure how could I use python to achieve it.

Thank you,  
Xinyue

---

<div class="post-metadata">

### Author: ![lassoan](https://discourse.itk.org/user_avatar/discourse.itk.org/lassoan/32/27_2.png) [@lassoan](https://discourse.itk.org/u/lassoan)
#### Post date: [August 16, 2023, 10:41am UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/7 "2023-08-16T10:41:14Z")

</div>

You can rigidly transform an image without resampling, just by updating the image geometry information (image origin, spacing, axis directions). This is what 3D Slicer does, too, when you harden a rigid transform on an image.

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [August 16, 2023, 1:27pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/8 "2023-08-16T13:27:45Z")

</div>

That can be accomplished in ITK via [TransformGeometryImageFilter](https://itk.org/Doxygen/html/classitk_1_1TransformGeometryImageFilter.html), see [here](https://discourse.itk.org/t/cuberille-mesh/6130/7).

---

<div class="post-metadata">

### Author: ![Xinyue](https://discourse.itk.org/user_avatar/discourse.itk.org/xinyue/32/3339_2.png) [@Xinyue](https://discourse.itk.org/u/Xinyue)
#### Post date: [August 16, 2023, 2:48pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/9 "2023-08-16T14:48:08Z")

</div>

Hi @lassoan and @dzenanz ,

Thank you for your reply!

The document mentioned that it transform the image “physically” by changing the metadata. May I ask that after GetArrayFromImage(), will this change keeps in the numpy array as well?

Thank you,  
Xinyue

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [August 16, 2023, 2:54pm UTC](https://discourse.itk.org/t/split-deformation-vector-field-into-two-deformation-vector-fields-and-apply-to-image-in-order/6148/10 "2023-08-16T14:54:54Z")

</div>

As numpy array throws away metadata, it will be unchanged by this filter.
