RTK/Examples/MCCBCTReconstruction: Difference between revisions

From Openrtk
Jump to navigation Jump to search
No edit summary
(13 intermediate revisions by 2 users not shown)
Line 1: Line 1:
RTK provides the necessary tools to reconstruct an image with motion compensation. The implementation is based on two articles that we have published ([http://www.creatis.insa-lyon.fr/site/fr/publications/RIT-09 article 1] and [http://www.creatis.insa-lyon.fr/site/fr/publications/RIT-09b article 2]) but only the FDK-based motion-compensated CBCT reconstruction (analytic algorithm in article 1) and without optimization (very slow reconstruction compared to article 2). You should read the articles to understand the basics of the algorithm before trying to use the software.
RTK provides the necessary tools to reconstruct an image with motion compensation. The implementation is based on two articles that we have published ([https://hal.archives-ouvertes.fr/hal-00443440 article 1] and [https://hal.archives-ouvertes.fr/hal-01967313 article 2]) but only the FDK-based motion-compensated CBCT reconstruction (analytic algorithm in article 1) and without optimization (very slow reconstruction compared to article 2). You should read the articles to understand the basics of the algorithm before trying to use the software.


The algorithm requires a set of projection images with the associated RTK geometry, the respiratory phase of each projection image and the 4D motion vector field over a respiratory cycle in the cone-beam coordinate system. Each piece of data is described in more details below. It is assumed that we have a breathing motion that is cyclic and similar to that described by the vector field. Note that you could modify the code and create your own motion model if you want to, in which case you should probably [http://www.openrtk.org/RTK/project/contactus.html contact us].
The algorithm requires a set of projection images with the associated RTK geometry, the respiratory phase of each projection image and the 4D motion vector field over a respiratory cycle in the cone-beam coordinate system. Each piece of data is described in more details below and can be downloaded using [https://data.kitware.com/#collection/5a7706878d777f0649e04776 Girder]. It is assumed that we have a breathing motion that is cyclic and similar to that described by the vector field. Note that you could modify the code and create your own motion model if you want to, in which case you should probably [http://www.openrtk.org/RTK/project/contactus.html contact us].


= Projection images =
= Projection images =


This example is illustrated with a set of projection images of the [http://www.creatis.insa-lyon.fr/rio/popi-model_original_page POPI patient]. This dataset has been used in the first previously-mentioned article. You can [http://midas3.kitware.com/midas/download/item/208648/projections.tgz download the projections] and the required tables of the Elekta database, [http://midas3.kitware.com/midas/download/item/208650/FRAME.DBF FRAME.DBF] and [http://midas3.kitware.com/midas/download/item/208649/IMAGE.DBF IMAGE.DBF], from MIDAS. The dataset is first used to reconstruct a blurry image:
This example is illustrated with a set of projection images of the [http://www.creatis.insa-lyon.fr/rio/popi-model_original_page POPI patient]. This dataset has been used in the first previously-mentioned article. You can [https://data.kitware.com/api/v1/item/5be99af88d777f2179a2e144/download download the projections] and the required tables of the Elekta database, [https://data.kitware.com/api/v1/item/5be99a068d777f2179a2cf4f/download FRAME.DBF] and [https://data.kitware.com/api/v1/item/5be99a078d777f2179a2cf65/download IMAGE.DBF]. The dataset is first used to reconstruct a blurry image:


<source lang="bash">
<source lang="bash">
Line 39: Line 39:
= Deformation vector field =
= Deformation vector field =


The next piece of data is a 4D deformation vector field that describes a respiratory cycle. Typically, it can be obtained from the 4D planning CT with deformable image registration. Here, I have used [http://elastix.bigr.nl/ Elastix] with the [http://elastix.bigr.nl/wiki/index.php/Par0016 new sliding module] developed by Vivien Delmon. The registration uses a [http://midas3.kitware.com/midas/download/item/208653/patient_50.mha patient mask] (red+green) and a [http://midas3.kitware.com/midas/download/item/208652/mm_50.mha motion mask] (red) as described in [http://www.creatis.insa-lyon.fr/site/fr/publications/VAND-12 Jef's publication]:
The next piece of data is a 4D deformation vector field that describes a respiratory cycle. Typically, it can be obtained from the 4D planning CT with deformable image registration. Here, I have used [http://elastix.bigr.nl/ Elastix] with the [http://elastix.bigr.nl/wiki/index.php/Par0016 sliding module] developed by Vivien Delmon. The registration uses a [https://data.kitware.com/api/v1/item/5be99a408d777f2179a2dde8/download patient mask] (red+green) and a [https://data.kitware.com/api/v1/item/5be99a088d777f2179a2cf6f/download motion mask] (red) as described in [http://www.creatis.insa-lyon.fr/site/fr/publications/VAND-12 Jef's publication]:


[[File:mm.jpg]]
[[File:mm.jpg]]




The registration can easily be scripted, here with bash:
The registration can easily be scripted, here with bash, where each phase image of the POPI 4D CT has been stored in files 00.mhd to 50.mhd:


<source lang="bash">
<source lang="bash">
Line 50: Line 50:
do
do
   mkdir $i
   mkdir $i
   elastix -f 50.mhd -m $i.mhd -out $i -labels mm_50.mha -fMask patient_50.mha -p Par0016.multibsplines.lung.sliding.txt
   elastix -f 50.mhd \
          -m $i.mhd \
          -out $i \
          -labels mm_50.mha \
          -fMask patient_50.mha \
          -p Par0016.multibsplines.lung.sliding.txt
done  
done  
</source>
</source>


Note that the reference frame is phase 50% and it is registered to each phase from 0% to 90%. One subtle step is that the vector field is a displacement vector field, i.e., each vector is the displacement of the point a the location of the arrow. Since I ran the registration on the 4D planning CT, the coordinate system is not that of the cone-beam CT. In order to produce the vector field in the cone-beam coordinate system, I have used the following bash script that combines transformix and several "clitk" tools that are provided along with VV:
Deformable Image Registration is a complex and long process so you will have to be patient here. Note that the reference frame is phase 50% and it is registered to each phase from 0% to 90%. One subtle step is that the vector field is a displacement vector field, i.e., each vector is the local displacement of the point at its location. Since I ran the registration on the 4D planning CT, the coordinate system is not that of the cone-beam CT. In order to produce the vector field in the cone-beam coordinate system, I have used the following bash script that combines transformix and several "clitk" tools that are provided along with VV:


<source lang="bash">
<source lang="bash">
# Create 4x4 matrix that describes the CT to CBCT change of coordinate system.
# This matrix is a combination of the knowledge of the isocenter position / axes orientation
# and a rigid alignment that has been performed with Elastix
echo "-0.0220916855767852  0.9996655273534405 -0.0134458487848415 -83.6625731437426197"  >CT_CBCT.mat
echo " 0.0150924269790251 -0.0131141301144939 -0.9998000991394341  -4.0763571826687057" >>CT_CBCT.mat
echo " 0.9996420239647088  0.0222901999207823  0.0147976657359281  77.8903364738220034" >>CT_CBCT.mat
echo " 0.0000000000000000  0.0000000000000000  0.0000000000000000  1.0000000000000000" >>CT_CBCT.mat
# Transform 4x4 matrix that describes the transformation
# Transform 4x4 matrix that describes the transformation
# from planning CT to CBCT to a vector field
# from planning CT to CBCT to a vector field
clitkMatrixTransformToVF --like 50.mhd \
clitkMatrixTransformToVF --like 50.mhd \
                         --input CT_CBCT.mat \
                         --matrix CT_CBCT.mat \
                         --output CT_CBCT.mha
                         --output CT_CBCT.mha


# Inverse transformation remove upper slices that are ouside the planning CT
# Inverse transformation. Also remove upper slices that are outside the
# CBCT_CT.mat is the inverse of CT_CBCT.mha
# planning CT CBCT_CT.mat is the inverse of CT_CBCT.mha
clitkMatrixInverse -i CT_CBCT.mat \
                  -o CBCT_CT.mat
clitkMatrixTransformToVF --origin -127.5,-107.5,-127.5 \
clitkMatrixTransformToVF --origin -127.5,-107.5,-127.5 \
                         --spacing 1,1,1 \
                         --spacing 1,1,1 \
                         --size 256,236,256 \
                         --size 256,236,256 \
                         --input CBCT_CT.mat \
                         --matrix CBCT_CT.mat \
                         --output CBCT_CT.mha
                         --output CBCT_CT.mha


Line 79: Line 94:
               -tp $i/TransformParameters.0.txt \
               -tp $i/TransformParameters.0.txt \
               -def all -threads 16
               -def all -threads 16
   clitkComposeVF -i CBCT_CT.mha \
   clitkComposeVF --input1 CBCT_CT.mha \
                 -j $i/deformationField.mhd \
                 --input2 $i/deformationField.mhd \
                 -o $i/deformationField.mhd
                 --output $i/deformationField.mhd
   clitkComposeVF -i $i/deformationField.mhd \
   clitkComposeVF --input1 $i/deformationField.mhd \
                 -j CT_CBCT.mha \
                 --input2 CT_CBCT.mha \
                 -o $i/deformationField.mhd
                 --output $i/deformationField.mhd
done
done
</source>
</source>


This is a bit complicated and there are probably other ways of doing this. For example, Vivien has resampled the planning CT frames on the CBCT coordinate system before doing the registrations, in which case you do not need to do all this. Just pick one of your choice but you should obtain at the end a 4D vector field that is nicely displayed on top of the fdk.mha file that RTK produce such as this (the vector field is downsampled and displayed with VV):
This is a bit complicated and there are probably other ways of doing this. For example, Vivien has resampled the planning CT frames on the CBCT coordinate system before doing the registrations, in which case you do not need to do all this. Just pick one of your choice but motion-compensated CBCT reconstruction requires a 4D vector field that is nicely displayed on top of a CBCT image, for example the fdk.mha that has been produced in the first step (the vector field is downsampled and displayed with VV):


[[File:vf.gif]]
[[File:vf.gif]]
The elastix output files and the transformed 4D DVF are available [https://data.kitware.com/api/v1/item/5be99a058d777f2179a2cf42/download here].


= Respiratory signal =
= Respiratory signal =
Line 97: Line 114:


<source lang="bash">
<source lang="bash">
rtkamsterdamshroud -p . -r '.*.his' -o shroud.mha -u 650
rtkamsterdamshroud --path . \
rtkextractshroudsignal -i shroud.mha -o signal.txt
                  --regexp '.*.his' \
                  --output shroud.mha \
                  --unsharp 650
rtkextractshroudsignal --input shroud.mha \
                      --output signal.txt
</source>
</source>


in combination with Matlab post-processing to get the phase signal. Note that the phase must go from 0 to 1 where 0.3 corresponds to 30% in the respiratory cycle, i.e., frame 3 if you have a 10-frames 4D DVF or frame 6 if you have a 20-frames 4D DVF. The [http://midas3.kitware.com/midas/download/item/208654/signal.txt resulting phase] is in green on top of the blue respiratory signal and the detected end-exhale peaks:
in combination with Matlab post-processing to get the phase signal. Note that the phase must go from 0 to 1 where 0.3 corresponds to 30% in the respiratory cycle, i.e., frame 3 if you have a 10-frames 4D DVF or frame 6 if you have a 20-frames 4D DVF. The [https://data.kitware.com/api/v1/item/5be99af98d777f2179a2e160/download resulting phase] is in green on top of the blue respiratory signal and the detected end-exhale peaks:


[[File:signal.gif]]
[[File:signal.gif]]
Line 107: Line 128:
= Motion-compensated cone-beam CT reconstruction =
= Motion-compensated cone-beam CT reconstruction =


We now have all the pieces to do a motion-compensated reconstruction. Gathering all these pieces is probably the main difficulty. If you have set everything correctly, you now only have to run the following command
We now have all the pieces to do a motion-compensated reconstruction. Gathering all these pieces is probably the main difficulty. If you have set everything correctly, you now only have to run the following commands
 
<source lang="bash">
# Reconstruct from all projection images with motion compensation
rtkfdk \
  -p . \
  -r .*.his \
  -o fdk.mha \
  -g geometry.rtk \
  --hann 0.5 \
  --pad 1.0 \
  --signal sphase.txt \
  --dvf deformationField_4D.mhd
 
# Keep only the field-of-view of the image
rtkfieldofview \
  --reconstruction fdk.mha \
  --output fdk.mha \
  --geometry geometry.rtk \
  --path . \
  --regexp '.*.his'
</source>
 
At this stage, you should have the time to drink a coffee (several for some computers). A neat solution to appreciate the improvement is to toggle between uncorrected and motion-compensated reconstruction:
 
[[File:blurred_vs_mc.gif]]
 
You should be aware that we have constructed the 4D vector field with phase 50% as a reference. Therefore, the image is compensated towards that position. However, any other phase can be reconstructed by modifying the reference image of your 4D vector field. We usually like to reconstruct towards the time-average position.

Revision as of 07:21, 30 April 2019

RTK provides the necessary tools to reconstruct an image with motion compensation. The implementation is based on two articles that we have published (article 1 and article 2) but only the FDK-based motion-compensated CBCT reconstruction (analytic algorithm in article 1) and without optimization (very slow reconstruction compared to article 2). You should read the articles to understand the basics of the algorithm before trying to use the software.

The algorithm requires a set of projection images with the associated RTK geometry, the respiratory phase of each projection image and the 4D motion vector field over a respiratory cycle in the cone-beam coordinate system. Each piece of data is described in more details below and can be downloaded using Girder. It is assumed that we have a breathing motion that is cyclic and similar to that described by the vector field. Note that you could modify the code and create your own motion model if you want to, in which case you should probably contact us.

Projection images

This example is illustrated with a set of projection images of the POPI patient. This dataset has been used in the first previously-mentioned article. You can download the projections and the required tables of the Elekta database, FRAME.DBF and IMAGE.DBF. The dataset is first used to reconstruct a blurry image:

# Convert Elekta database to RTK geometry
rtkelektasynergygeometry \
  -o geometry.rtk \
  -f FRAME.DBF \
  -i IMAGE.DBF \
  -u 1.3.46.423632.141000.1169042526.68

# Reconstruct from all projection images without any motion compensation
rtkfdk \
  -p . \
  -r .*.his \
  -o fdk.mha \
  -g geometry.rtk \
  --hann 0.5 \
  --pad 1.0

# Keep only the field-of-view of the image
rtkfieldofview \
  --reconstruction fdk.mha \
  --output fdk.mha \
  --geometry geometry.rtk \
  --path . \
  --regexp '.*.his'

You should obtain something like that with VV:

Blurred.jpg

Deformation vector field

The next piece of data is a 4D deformation vector field that describes a respiratory cycle. Typically, it can be obtained from the 4D planning CT with deformable image registration. Here, I have used Elastix with the sliding module developed by Vivien Delmon. The registration uses a patient mask (red+green) and a motion mask (red) as described in Jef's publication:

Mm.jpg


The registration can easily be scripted, here with bash, where each phase image of the POPI 4D CT has been stored in files 00.mhd to 50.mhd:

for i in $(seq -w 0 10 90)
do
  mkdir $i
  elastix -f 50.mhd \
          -m $i.mhd \
          -out $i \
          -labels mm_50.mha \
          -fMask patient_50.mha \
          -p Par0016.multibsplines.lung.sliding.txt
done

Deformable Image Registration is a complex and long process so you will have to be patient here. Note that the reference frame is phase 50% and it is registered to each phase from 0% to 90%. One subtle step is that the vector field is a displacement vector field, i.e., each vector is the local displacement of the point at its location. Since I ran the registration on the 4D planning CT, the coordinate system is not that of the cone-beam CT. In order to produce the vector field in the cone-beam coordinate system, I have used the following bash script that combines transformix and several "clitk" tools that are provided along with VV:

# Create 4x4 matrix that describes the CT to CBCT change of coordinate system.
# This matrix is a combination of the knowledge of the isocenter position / axes orientation
# and a rigid alignment that has been performed with Elastix
echo "-0.0220916855767852  0.9996655273534405 -0.0134458487848415 -83.6625731437426197"  >CT_CBCT.mat
echo " 0.0150924269790251 -0.0131141301144939 -0.9998000991394341  -4.0763571826687057" >>CT_CBCT.mat
echo " 0.9996420239647088  0.0222901999207823  0.0147976657359281  77.8903364738220034" >>CT_CBCT.mat
echo " 0.0000000000000000  0.0000000000000000  0.0000000000000000   1.0000000000000000" >>CT_CBCT.mat

# Transform 4x4 matrix that describes the transformation
# from planning CT to CBCT to a vector field
clitkMatrixTransformToVF --like 50.mhd \
                         --matrix CT_CBCT.mat \
                         --output CT_CBCT.mha

# Inverse transformation. Also remove upper slices that are outside the
# planning CT CBCT_CT.mat is the inverse of CT_CBCT.mha
clitkMatrixInverse -i CT_CBCT.mat \
                   -o CBCT_CT.mat
clitkMatrixTransformToVF --origin -127.5,-107.5,-127.5 \
                         --spacing 1,1,1 \
                         --size 256,236,256 \
                         --matrix CBCT_CT.mat \
                         --output CBCT_CT.mha

# Go over each elastix output file, generate the vector field with
# transformix and compose with the two rigid vector fields
for i in $(seq -w 0 10 90)
do
  transformix -in 50.mhd \
              -out $i \
              -tp $i/TransformParameters.0.txt \
              -def all -threads 16
  clitkComposeVF --input1 CBCT_CT.mha \
                 --input2 $i/deformationField.mhd \
                 --output $i/deformationField.mhd
  clitkComposeVF --input1 $i/deformationField.mhd \
                 --input2 CT_CBCT.mha \
                 --output $i/deformationField.mhd
done

This is a bit complicated and there are probably other ways of doing this. For example, Vivien has resampled the planning CT frames on the CBCT coordinate system before doing the registrations, in which case you do not need to do all this. Just pick one of your choice but motion-compensated CBCT reconstruction requires a 4D vector field that is nicely displayed on top of a CBCT image, for example the fdk.mha that has been produced in the first step (the vector field is downsampled and displayed with VV):

Vf.gif

The elastix output files and the transformed 4D DVF are available here.

Respiratory signal

The motion model requires that we associate each projection image with one frame of the 4D vector field. We used the Amsterdam shroud solution of Lambert Zijp (described here) which is implemented in RTK

rtkamsterdamshroud --path . \
                   --regexp '.*.his' \
                   --output shroud.mha \
                   --unsharp 650
rtkextractshroudsignal --input shroud.mha \
                       --output signal.txt

in combination with Matlab post-processing to get the phase signal. Note that the phase must go from 0 to 1 where 0.3 corresponds to 30% in the respiratory cycle, i.e., frame 3 if you have a 10-frames 4D DVF or frame 6 if you have a 20-frames 4D DVF. The resulting phase is in green on top of the blue respiratory signal and the detected end-exhale peaks:

Signal.gif

Motion-compensated cone-beam CT reconstruction

We now have all the pieces to do a motion-compensated reconstruction. Gathering all these pieces is probably the main difficulty. If you have set everything correctly, you now only have to run the following commands

# Reconstruct from all projection images with motion compensation
rtkfdk \
  -p . \
  -r .*.his \
  -o fdk.mha \
  -g geometry.rtk \
  --hann 0.5 \
  --pad 1.0 \
  --signal sphase.txt \
  --dvf deformationField_4D.mhd

# Keep only the field-of-view of the image
rtkfieldofview \
  --reconstruction fdk.mha \
  --output fdk.mha \
  --geometry geometry.rtk \
  --path . \
  --regexp '.*.his'

At this stage, you should have the time to drink a coffee (several for some computers). A neat solution to appreciate the improvement is to toggle between uncorrected and motion-compensated reconstruction:

Blurred vs mc.gif

You should be aware that we have constructed the 4D vector field with phase 50% as a reference. Therefore, the image is compensated towards that position. However, any other phase can be reconstructed by modifying the reference image of your 4D vector field. We usually like to reconstruct towards the time-average position.