RTK/Examples/4DROOSTERReconstruction

RTK provides a tool to reconstruct a 3D + time image which does not require explicit motion estimation. Instead, it uses total variation regularization along both space and time. The implementation is based on a paper that we have published (article). You should read the article 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 a motion mask in which the region of the space where movement is expected is set to 1 and the rest is set to 0. Each piece of data is described in more details below and can be downloaded from MIDAS. It is assumed that the breathing motion is periodic, which implies that the mechanical state of the chest depends only on the respiratory phase.

Projection images

This example is illustrated with a set of projection images of the POPI patient. 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 a 3D volume from all projection images. It implicitly assumes that the patient's chest is static, which is wrong. Therefore the reconstruction will be blurry
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:

Motion mask

The next piece of data is a 3D motion mask: a volume that contains only zeros, where no movement is expected to occur, and ones, where movement is expected. Typically, during breathing, the whole chest moves, so it is safe to use the patient mask (red+green). However, restricting the movement to a smaller region, e.g. the rib cage, can help reconstruct this region more accurately and mitigate artifacts, so we might want to use the rib cage mask (red):

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:

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:

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.