[3]:
# -*- coding: utf-8 -*-
# Copyright 2023 United Kingdom Research and Innovation
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authored by: Tommi Heikkilä (University of Helsinki, Finland),
# Edited by: Margaret Duff (STFC-UKRI)
TV-regularized reconstruction of the dynamix STEMPO dataset in CIL#
Based mostly on the work carried out during the CIL Hackathon in Cambridge, UK 20-23 March 2023.
Created by Tommi Heikkilä, University of Helsinki, Finland.
Special thanks to the CIL team and Christoph Kolbitsch (PTB) for all the help.
Computations performed using CIL 23.0.1.
General outline#
Setup: CIL, load the data, the geometry etc.
Start with static data and single time step reconstruction
Move to dynamic data and single time step reconstruction
Split dynamic data into multiple time steps (batches)
Reconstruct each at the same time using FISTA and (spatial) TV regularization
(Inefficient temporal + spatial TV implementation)
Future goal#
Temporal (anisotropic) TV regularization (i.e regularize across time steps) for controlled temporal evolution. Preferably in a more efficient manner, similar to Papoutsellis E et al. (2021) Core Imaging Library - Part II: multichannel reconstruction for dynamic and spectral tomography. Phil. Trans. R. Soc. A 20200193.
Note: it is very likely that the approach taken here (Partitioner) is not suitable for this.
License#
This code is licensed under the Apache License, Version 2.0 which is availalbe at: http://www.apache.org/licenses/LICENSE-2.0
The STEMPO Dataset is licensed under Creative Commons Attributiom 4.0 International: https://creativecommons.org/licenses/by/4.0/legalcode
[4]:
# cil imports
from cil.framework import ImageData, ImageGeometry
from cil.framework import AcquisitionGeometry, AcquisitionData, BlockDataContainer
from cil.processors import Slicer, AbsorptionTransmissionConverter, TransmissionAbsorptionConverter
from cil.optimisation.functions import IndicatorBox, LeastSquares, ZeroFunction
from cil.optimisation.algorithms import CGLS, SIRT, FISTA
from cil.optimisation.operators import BlockOperator, GradientOperator, IdentityOperator, FiniteDifferenceOperator
from cil.plugins.astra.operators import ProjectionOperator
from cil.plugins.astra.processors import FBP
from cil.recon import FDK
from cil.plugins import TomoPhantom
from cil.utilities import dataexample
from cil.utilities.display import show2D, show1D, show_geometry
# External imports
import numpy as np
import matplotlib.pyplot as plt
import logging
import scipy.io
import os
import mat73 # This makes loading MATLAB v7.3 files much easier # can be installed via pip
import cil
print(cil.__version__)
23.0.1
[5]:
# Set logging level for CIL processors:
logging.basicConfig(level=logging.WARNING)
cil_log_level = logging.getLogger('cil.processors')
cil_log_level.setLevel(logging.INFO)
1. Load STEMPO data#
Data is available in Zenodo: https://zenodo.org/record/7147139 / https://doi.org/10.5281/zenodo.7147139
Let’s start with the static data to make sure the geometry is set up correctly
Note: it may be worth using smaller resolution data (larger binning value) than \(b=4\).
[6]:
filename = "stempo_static_2d_b4.mat"
data = mat73.loadmat(filename)
CtData = data["CtData"]
sinogram = np.float32(CtData["sinogram"])
# Load parameters
param = CtData["parameters"]
angles = np.float32(param["angles"])
[7]:
print(f'Sinogram size is {sinogram.shape}')
Sinogram size is (360, 560)
1.1 Set up one large geometry#
[8]:
# Acquisition geometry parameters
binningFactor = param["binningPost"] # Data is binned
SOD = param["distanceSourceOrigin"]
print(f'SOD: {SOD:.3f} mm')
SDD = param["distanceSourceDetector"]
print(f'SDD: {SDD:.3f} mm')
ODD = SDD - SOD # Origin-detector-distance
print(f'ODD: {ODD:.3f} mm')
pixelSize = param["pixelSizePost"] # Detector pixel width after binning
print(f'Detector pixel size: {pixelSize:.3f} mm')
SOD: 410.660 mm
SDD: 553.740 mm
ODD: 143.080 mm
Detector pixel size: 0.200 mm
Note that we can not change the rotational axis direction in fan-beam geometry so we simply change the sign of the angles
[9]:
ag = AcquisitionGeometry.create_Cone2D([0,-SOD], [0,ODD], units='mm',detector_direction_x=[1, 0])
num_angles, num_pixels = sinogram.shape
ag.set_panel(num_pixels, pixel_size=pixelSize, origin='bottom-left')
ag.set_angles(-angles)
ag.set_labels(('angle','horizontal'))
print(ag)
2D Cone-beam tomography
System configuration:
Source position: [ 0. , -410.66]
Rotation axis position: [0., 0.]
Detector position: [ 0. , 143.08]
Detector direction x: [1., 0.]
Panel configuration:
Number of pixels: [560 1]
Pixel size: [0.2 0.2]
Pixel origin: bottom-left
Channel configuration:
Number of channels: 1
Acquisition description:
Number of positions: 360
Angles 0-20 in degrees:
[ -0., -1., -2., -3., -4., -5., -6., -7., -8., -9., -10., -11.,
-12., -13., -14., -15., -16., -17., -18., -19.]
Distances in units: mm
[10]:
show_geometry(ag)
[10]:
<cil.utilities.display.show_geometry at 0x7fcc07ee6c20>
[11]:
# Setup image geometry
n = num_pixels
ig = ag.get_ImageGeometry()
ig.voxel_num_x = n
ig.voxel_num_y = n
print(ig)
Number of channels: 1
channel_spacing: 1.0
voxel_num : x560,y560
voxel_size : x0.14832231733304438,y0.14832231733304438
center : x0,y0
[12]:
# Setup data
stempoData = AcquisitionData(sinogram, deep_copy=False, geometry=ag)
[13]:
show2D(stempoData,size=(7,10))
[13]:
<cil.utilities.display.show2D at 0x7fcc07fc5900>
2. FBP of the static data#
[14]:
fbp = FBP(ig, ag, 'gpu')
[15]:
# run on the AcquisitionData
recon_fbp_static = fbp(stempoData)
[16]:
show2D(recon_fbp_static,size=(10,10))
[16]:
<cil.utilities.display.show2D at 0x7fcbffb16b30>
3. Repeat the previous steps with dynamic data#
Note: variable names are reused
[17]:
# Load data
filename = "stempo_seq8x45_2d_b4.mat"
data = mat73.loadmat(filename)
CtData = data["CtData"]
sinogram = np.float32(CtData["sinogram"])
# Load parameters
param = CtData["parameters"]
angles = np.float32(param["angles"])
binningFactor = param["binningPost"]
SOD = param["distanceSourceOrigin"]
print(f'SOD: {SOD:.3f} mm')
SDD = param["distanceSourceDetector"]
print(f'SDD: {SDD:.3f} mm')
ODD = SDD - SOD # Origin-detector-distance
print(f'ODD: {ODD:.3f} mm')
pixelSize = param["pixelSizePost"]
print(f'Detector pixel size: {pixelSize:.3f} mm')
# Acquisition geometry
ag = AcquisitionGeometry.create_Cone2D([0,-SOD], [0,ODD], units='mm',detector_direction_x=[1, 0])
num_angles, num_pixels = sinogram.shape
ag.set_panel(num_pixels, pixel_size=pixelSize, origin='bottom-left')
ag.set_angles(-angles)
ag.set_labels(('angle','horizontal'))
print(ag)
# Setup image geometry
n = 540
ig = ag.get_ImageGeometry()
ig.voxel_num_x = n
ig.voxel_num_y = n
print(ig)
SOD: 410.660 mm
SDD: 553.740 mm
ODD: 143.080 mm
Detector pixel size: 0.200 mm
2D Cone-beam tomography
System configuration:
Source position: [ 0. , -410.66]
Rotation axis position: [0., 0.]
Detector position: [ 0. , 143.08]
Detector direction x: [1., 0.]
Panel configuration:
Number of pixels: [560 1]
Pixel size: [0.2 0.2]
Pixel origin: bottom-left
Channel configuration:
Number of channels: 1
Acquisition description:
Number of positions: 360
Angles 0-20 in degrees:
[ -0., -8., -16., -24., -32., -40., -48., -56., -64., -72.,
-80., -88., -96., -104., -112., -120., -128., -136., -144., -152.]
Distances in units: mm
Number of channels: 1
channel_spacing: 1.0
voxel_num : x540,y540
voxel_size : x0.14832231733304438,y0.14832231733304438
center : x0,y0
[18]:
# Setup data
stempoData = AcquisitionData(sinogram, deep_copy=False, geometry=ag)
show2D(stempoData, size=(10,10))
[18]:
<cil.utilities.display.show2D at 0x7fcbffb98e50>
The sinogram looks rather busy since it contains (sparse set of) projections from 8 full rotations around the object.
[19]:
plt.imshow(sinogram, cmap='gray', aspect=3)
[19]:
<matplotlib.image.AxesImage at 0x7fcbffbcb850>
3.1 Witness the movement artefacts#
Using the adjoint of the projection operator we can see some of the movement artefacts:
[20]:
# Create projection operator using Astra-Toolbox.
A = ProjectionOperator(ig, ag, 'gpu')
bp = A.adjoint(stempoData) # Backgprojection
show2D([ bp], title=[ 'Adjoint reconstruction'], size=(10,10), num_cols=2)
[20]:
<cil.utilities.display.show2D at 0x7fcbf519dcc0>
3.2 FBP reconstruction of a dynamic data#
Compare the back projection from the dynamic data from the previous reconstruction from static data:
[21]:
fbp = FBP(ig, ag, 'gpu')
# run on the AcquisitionData
recon_fbp = fbp(stempoData)
show2D([recon_fbp_static, recon_fbp], title=['Static reconstruction', 'Dynamic reconstruction'], size=(10,10), num_cols=2)
[21]:
<cil.utilities.display.show2D at 0x7fcbf4f68370>
That looks terrible. Let’s remove negative pixel values by fixing the range of the plot:
[22]:
show2D(recon_fbp, size=(10,10), fix_range=(0,0.05))
[22]:
<cil.utilities.display.show2D at 0x7fcc08099ed0>
Now it only looks bad. We can see notable motion artefacts in the middle and a lot of sparse angle artefacts since we only have 45 unique projection directions.
3.3 Iterative methods have the same artefacts#
First use CGLS to solve the optimisation problem
\[\min_x\|Ax-y\|_2^2\]
where \(x\) is the image to reconstruct, \(A\) is the projection operator, and \(y\) is the observed sinogram.
[23]:
# initialisation of the optimisation algorithm - zero array in this case
initial = ig.allocate(0)
# setup CGLS
cgls = CGLS(initial=initial,
operator=A,
data=stempoData,
max_iteration = 100,
update_objective_interval = 10 )
Running for 50 iterations gives a similar reconstruction to the FBP:
[24]:
cgls.run(50,verbose=True)
# get and visualise the results
recon_cgls = cgls.solution
Iter Max Iter Time/Iter Objective
[s]
0 100 0.000 6.73189e+04
10 100 0.061 3.45551e+03
20 100 0.060 3.43862e+03
30 100 0.059 3.43665e+03
40 100 0.068 3.43600e+03
50 100 0.080 3.43566e+03
-------------------------------------------------------
50 100 0.080 3.43566e+03
Stop criterion has been reached.
We can also try running for a smaller number of iterations to regularise by early stopping:
[25]:
# initialisation of the optimisation algorithm - zero array in this case
initial = ig.allocate(0)
# setup CGLS
cgls = CGLS(initial=initial,
operator=A,
data=stempoData,
max_iteration = 100,
update_objective_interval = 10 )
cgls.run(5,verbose=True)
# get and visualise the results
recon_cgls_early_stopping = cgls.solution
Iter Max Iter Time/Iter Objective
[s]
0 100 0.000 6.73189e+04
-------------------------------------------------------
5 100 0.117
Stop criterion has been reached.
[26]:
show2D([recon_fbp, recon_cgls, recon_cgls_early_stopping], ['FBP', 'CGLS', 'CGLS_early_stopping'], \
num_cols=3, size=(10,10), origin='upper-left', fix_range=(0,0.05))
[26]:
<cil.utilities.display.show2D at 0x7fcbf4c28d30>
3.4 Consider only one time step (out of 8)#
[27]:
# Limited angles
# We cut a small section out of the sinogram and angles
numAngles = 45
anglesSparse = angles[0:numAngles]
sinoSparse = sinogram[0:numAngles,:]
print(f'Sinogram size is now {sinoSparse.shape}')
agSparse = ag.copy()
agSparse.set_angles(-anglesSparse)
stempoSparse = AcquisitionData(sinoSparse, deep_copy=False, geometry=agSparse)
show2D(stempoSparse, size=(10,10))
Sinogram size is now (45, 560)
[27]:
<cil.utilities.display.show2D at 0x7fcbf4c41060>
[28]:
fbp = FBP(ig, agSparse, 'gpu')
# run on the AcquisitionData
sparse_fbp = fbp(stempoSparse)
show2D(sparse_fbp, size=(10,10), fix_range=(0,0.05), title='t=1, 45 angles FBP')
[28]:
<cil.utilities.display.show2D at 0x7fcbf4f7a0e0>
Some movement artefacts are still visible but the situations has improved. Let’s try iterative algorithm (cgls) and consider regularisation by early stopping:
[29]:
A = ProjectionOperator(ig, agSparse, 'gpu')
# initial estimate - zero array in this case
initial = ig.allocate(0)
# setup CGLS
cgls = CGLS(initial=initial,
operator=A,
data=stempoSparse,
max_iteration = 100,
update_objective_interval = 10 )
cgls.run(50,verbose=True)
# get and visualise the results
sparse_cgls = cgls.solution
# initialisation of the optimisation algorithm - zero array in this case
initial = ig.allocate(0)
# setup CGLS
cgls = CGLS(initial=initial,
operator=A,
data=stempoSparse,
max_iteration = 100,
update_objective_interval = 10 )
cgls.run(5,verbose=True)
# get and visualise the results
sparse_cgls_early_stopping = cgls.solution
show2D([sparse_fbp, sparse_cgls, sparse_cgls_early_stopping], ['FBP', 'CGLS', 'CGLS_early_stopping'], \
num_cols=3, size=(10,10), origin='upper-left', fix_range=(0,0.05))
Iter Max Iter Time/Iter Objective
[s]
0 100 0.000 8.42717e+03
10 100 0.078 1.65551e+01
20 100 0.070 1.36376e+01
30 100 0.069 1.30771e+01
40 100 0.070 1.28557e+01
50 100 0.062 1.27849e+01
-------------------------------------------------------
50 100 0.062 1.27849e+01
Stop criterion has been reached.
Iter Max Iter Time/Iter Objective
[s]
0 100 0.000 8.42717e+03
-------------------------------------------------------
5 100 0.026
Stop criterion has been reached.
[29]:
<cil.utilities.display.show2D at 0x7fcbf4f680d0>
Try Tikhonov regularisation:
\[\min_x\|Ax-y\|_2^2+\alpha\|x\|_2^2\]
where \(\alpha\) is a regularisation parameter that balances the data discrepency term and the regularisation term:
[30]:
# Tikhonov regularization with one time step
L = IdentityOperator(ig)
alpha = 0.7
operator_block = BlockOperator(A, alpha*L)
zero_data = L.range.allocate(0)
data_block = BlockDataContainer(stempoSparse, zero_data)
#setup CGLS with the Block Operator and Block DataContainer
initial = ig.allocate(0)
cgls_tikh = CGLS(initial=initial, operator=operator_block, data=data_block, update_objective_interval = 10, tolerance=1e-05)
cgls_tikh.max_iteration = 1000
#run the algorithm until convergence
cgls_tikh.run(1000)
show2D(cgls_tikh.solution, size=(10,10), fix_range=(0,0.15))
Iter Max Iter Time/Iter Objective
[s]
0 1000 0.000 8.42717e+03
10 1000 0.059 3.78615e+01
20 1000 0.075 3.55106e+01
30 1000 0.081 3.52775e+01
40 1000 0.083 3.52379e+01
50 1000 0.086 3.52269e+01
60 1000 0.088 3.52224e+01
70 1000 0.087 3.52186e+01
80 1000 0.087 3.52170e+01
90 1000 0.086 3.52167e+01
100 1000 0.085 3.52175e+01
110 1000 0.084 3.52195e+01
120 1000 0.084 3.52204e+01
130 1000 0.085 3.52230e+01
140 1000 0.083 3.52247e+01
150 1000 0.080 3.52255e+01
160 1000 0.077 3.52256e+01
170 1000 0.075 3.52266e+01
180 1000 0.073 3.52284e+01
190 1000 0.071 3.52284e+01
200 1000 0.069 3.52284e+01
210 1000 0.067 3.52282e+01
220 1000 0.066 3.52285e+01
230 1000 0.066 3.52282e+01
240 1000 0.067 3.52266e+01
250 1000 0.068 3.52250e+01
260 1000 0.069 3.52230e+01
270 1000 0.070 3.52216e+01
280 1000 0.071 3.52198e+01
290 1000 0.072 3.52187e+01
300 1000 0.073 3.52180e+01
310 1000 0.073 3.52176e+01
320 1000 0.074 3.52175e+01
330 1000 0.074 3.52170e+01
340 1000 0.075 3.52169e+01
350 1000 0.075 3.52168e+01
360 1000 0.074 3.52166e+01
370 1000 0.073 3.52164e+01
380 1000 0.072 3.52164e+01
390 1000 0.071 3.52164e+01
400 1000 0.070 3.52162e+01
401 1000 0.070
Tolerance is reached: 1e-05
-------------------------------------------------------
401 1000 0.070
Stop criterion has been reached.
[30]:
<cil.utilities.display.show2D at 0x7fcbf509c190>
3.5 Just in case, test FISTA with one time step and non negativity constraint#
[31]:
A = ProjectionOperator(ig, agSparse, 'gpu')
LS = LeastSquares(A=A, b=stempoSparse)
G = IndicatorBox(lower=0.0)
fista_LS = FISTA(initial=ig.allocate(0), f=LS, g=G, max_iteration=1000, update_objective_interval=10)
[32]:
fista_LS.run(100)
Iter Max Iter Time/Iter Objective
[s]
0 1000 0.000 8.42717e+03
10 1000 0.021 7.46408e+01
20 1000 0.021 1.99271e+01
30 1000 0.021 1.60728e+01
40 1000 0.021 1.52739e+01
50 1000 0.021 1.49438e+01
60 1000 0.021 1.48031e+01
70 1000 0.021 1.47198e+01
80 1000 0.021 1.46630e+01
90 1000 0.023 1.46259e+01
100 1000 0.023 1.45980e+01
-------------------------------------------------------
100 1000 0.023 1.45980e+01
Stop criterion has been reached.
[33]:
LS_reco = fista_LS.solution
show2D(LS_reco, title="LS reconstruction", fix_range=(0,0.05), size=(10,10))
[33]:
<cil.utilities.display.show2D at 0x7fcbb7488b20>
4. Split into separate time steps#
We want to treat data mismatch of each time step separately but consider only one 2d + time object \(f\), i.e. something like
\[\Big\| \sum_{t=0}^{T-1} A f_t - m_t \Big\|^2 + \alpha R(f)\]
[34]:
## Split time steps
T = 8
from cil.framework.framework import Partitioner, DataContainer
from cil.optimisation.functions import L2NormSquared, L1Norm
from cil.optimisation.functions import MixedL21Norm
from cil.optimisation.functions import BlockFunction
from cil.optimisation.operators import ZeroOperator
stempoSplit = stempoData.partition(T, Partitioner.SEQUENTIAL)
[35]:
show2D(stempoSplit)
[35]:
<cil.utilities.display.show2D at 0x7fcbb7488970>
4.1 Build block (diagonal) operators#
Regular block operators and functions are block-like either in row or column direction, what we want is something like
\[\begin{split}A_{big} =
\left(\begin{array}{cccc}
A_0 & 0 & \dots & 0 \\
0 & A_1 & \dots & 0 \\
\vdots & & \ddots & \vdots \\
0 & 0 & \dots & A_{T-1}
\end{array}\right).\end{split}\]
Unfortunately at the moment CIL can not create such operators in efficient (sparse) way. So we have to do it inefficiently.
[36]:
projOp = ProjectionOperator(ig, stempoSplit.geometry)
projOpList = projOp.get_as_list() # List of the operators we want
# We need a zero-operator for padding (If the sizes changed for different time steps this would be pain)
# (Luckily they don't)
ZOp = ZeroOperator(projOpList[0].domain_geometry(), projOpList[0].range_geometry())
Alist = [ZOp for i in range(T*T)] # Long list of zero-operators
# Fill diagonal with proper operators
for t in range(T):
Alist[t*T+t] = projOpList[t]
# Obtain block-diagonal square operator
A = BlockOperator(*Alist, shape=(T,T))
bList = [stempoSplit.get_item(i) for i in range(T)]
b = BlockDataContainer(*bList)
[37]:
# Initial value using A
x0 = A.adjoint(b)
# Sanity check: we should see 8 different images
show2D(x0)
[37]:
<cil.utilities.display.show2D at 0x7fcbb798af20>
[38]:
# We initialize with zeros, which we have to do by iterating each image
for t in range(T):
x0[t].fill(0.0)
5. Iterative algorithms with multiple time steps#
5.1 FISTA with non-negativity constraint for each time step#
[39]:
# Nonnegativity constraint block function to act on each image separately
IBList = [IndicatorBox(lower=0.0) for i in range(T)]
IB = BlockFunction(*IBList)
# One large least squares term
F_LS = LeastSquares(A=A, b=stempoSplit)
# Set up FISTA algorithm
fista = FISTA(initial=x0, f=F_LS, g=IB, max_iteration=200, update_objective_interval=10)
[40]:
fista.run(100)
Iter Max Iter Time/Iter Objective
[s]
0 200 0.000 6.73190e+04
10 200 0.215 5.86625e+03
20 200 0.217 2.03110e+03
30 200 0.217 8.11439e+02
40 200 0.219 4.63374e+02
50 200 0.219 3.24651e+02
60 200 0.219 2.71682e+02
70 200 0.218 2.53289e+02
80 200 0.218 2.43930e+02
90 200 0.219 2.38423e+02
100 200 0.218 2.35406e+02
-------------------------------------------------------
100 200 0.218 2.35406e+02
Stop criterion has been reached.
[41]:
FISTArecn = fista.solution
show2D(FISTArecn, title=[f'FISTA, t={t+1}' for t in range(T)])
[41]:
<cil.utilities.display.show2D at 0x7fcbb675ddb0>
5.2 Let’s do some (spatial) TV regularization#
[42]:
from cil.plugins.ccpi_regularisation.functions import FGP_TV
[43]:
alpha = 0.05
TVlist = [FGP_TV(alpha=alpha, nonnegativity=True, device='gpu', isotropic=False) for t in range(T)]
TV = BlockFunction(*TVlist) # Unroll the list as before
fista_TV = FISTA(initial=x0, f=F_LS, g=TV, max_iteration=500, update_objective_interval=10)
[44]:
fista_TV.run(100)
Iter Max Iter Time/Iter Objective
[s]
0 500 0.000 6.73190e+04
10 500 0.532 5.88299e+03
20 500 0.423 2.07038e+03
30 500 0.390 8.61867e+02
40 500 0.371 5.19534e+02
50 500 0.361 3.84051e+02
60 500 0.354 3.32823e+02
70 500 0.348 3.15552e+02
80 500 0.345 3.06732e+02
90 500 0.342 3.01659e+02
100 500 0.339 2.98781e+02
-------------------------------------------------------
100 500 0.339 2.98781e+02
Stop criterion has been reached.
[45]:
FISTATVrecn = fista_TV.solution
show2D(FISTATVrecn, title=[f'FISTA+TV, t={t+1}' for t in range(T)])
[45]:
<cil.utilities.display.show2D at 0x7fcbb6f356c0>
This doesn’t seem too bad. However our block-diagonal approach does have its limits. Let’s push them!
6. Hacky temporal TV#
i.e. so many block diagonal operators…
For temporal TV we want to subtract subsequent time steps: \(f_t - f_{t+1}\). We can do this with a linear operator:
\[\begin{split}D_{t} =
\left(\begin{array}{ccccc}
I & -I & 0 & \dots & 0 \\
0 & I & -I & \dots & 0 \\
\vdots & & & \ddots & \vdots \\
0 & & \dots & I & -I \\
0 & & \dots & & 0
\end{array}\right)\end{split}\]
and then do $R(f) = :nbsphinx-math:`alpha `| D_t f |_2^2 $.
Problem is, we have very limited options for creating the operator such that the optimization methods still work. Luckily there is a way…
6.1 Let’s use MASSIVE block operators#
(This will be slow)
\[\begin{split}M =
\left(\begin{array}{cc}
A_{big} \\
\alpha D_t
\end{array}\right) =
\left(\begin{array}{ccc}
A_0 & \dots & 0 \\
& \ddots & \\
0 & \dots & A_{T-1} \\
\alpha I & -\alpha I & \\
& \ddots & \\
0 & \dots & 0
\end{array}\right)
, \quad B = \left(\begin{array}{cc}
b_0 \\
\vdots \\
b_{T-1} \\
0 \\
\vdots \\
0
\end{array}\right)\end{split}\]
and do least squares term with that
\[\| M f - B \|_2^2 = \| A_{big} f - b \|_2^2 + \alpha^2 \| D_t f \|_2^2.\]
This way our functional will be
\[\mathcal{F}(f) = \sum_{t=0}^{T-1} \left( \| A f_t - b_t \|_2^2 + \beta \, TV(f_t) \right) + \alpha^2 \| D_t f \|_2^2 .\]
[46]:
from cil.optimisation.operators import IdentityOperator
projOp = ProjectionOperator(ig, stempoSplit.geometry)
projOpList = projOp.get_as_list() # List of the operators we want
# We want the zero operator to have same domain and range as the forward operator projOp
ZOp1 = ZeroOperator(projOpList[0].domain_geometry(), projOpList[0].range_geometry())
Alist = [ZOp1 for i in range(T*T)] # Long list of zero-operators
# Fill diagonal with proper operators
for t in range(T):
Alist[t*T+t] = projOpList[t]
# Now we want to make each block operate from image domain to image domain!
ZOp2 = ZeroOperator(projOpList[0].domain_geometry(), projOpList[0].domain_geometry())
IOp = IdentityOperator(projOpList[0].domain_geometry(), projOpList[0].domain_geometry())
alpha = 5 # Regularization parameter for temporal TV
DtList = [ZOp2 for i in range(T*T)] # Long list of zero-operators
# Fill diagonal and off-diagonal with proper operators (last row is zeroes for simplicity)
for t in range(T-1):
DtList[t*T+t] = alpha*IOp
DtList[np.min((t,T-2))*T+t+1] = -alpha*IOp
# Use existing lists of operators and concanate
Mlist = [*Alist, *DtList]
# Obtain block-diagonal rectangular operator
M = BlockOperator(*Mlist, shape=(2*T,T))
bList = [stempoSplit.get_item(i) for i in range(T)]
x0List = [x0[i] for i in range(T)]
B = BlockDataContainer(*[*bList, *x0List]) # Concanate lists
MASSIVE_LS = LeastSquares(A=M, b=B)
[47]:
# Let's add spatial TV to the mix (luckily it can also include non-negativity)
beta = 0.05 # Regularization parameter for spatial TV
TV = BlockFunction(*[FGP_TV(alpha=beta, nonnegativity=True, device='gpu', isotropic=False) for t in range(T)])
[48]:
# We can set function g to be just non-negativity constraint (g=IB) or non-negativity + spatial TV (g=TV)
# Note that including spatial TV will make the algorithm even slower!
fista_MASSIVE = FISTA(initial=x0, f=MASSIVE_LS, g=TV, max_iteration=100, update_objective_interval=10)
So far so good. Fingers crossed it runs…
[49]:
fista_MASSIVE.run(50)
Iter Max Iter Time/Iter Objective
[s]
0 100 0.000 6.73190e+04
10 100 0.347 6.32259e+03
20 100 0.353 2.61179e+03
30 100 0.353 1.43551e+03
40 100 0.422 1.06185e+03
50 100 0.439 9.13366e+02
-------------------------------------------------------
50 100 0.439 9.13366e+02
Stop criterion has been reached.
[50]:
MASSIVErecn = fista_MASSIVE.solution
show2D(MASSIVErecn, title=[f'FISTA+TV, t={t+1}' for t in range(T)])
[50]:
<cil.utilities.display.show2D at 0x7fcbb5a38f10>
Comments#
It is possible to obtain (isotropic) temporal TV this way. Whether you actually should do it is another question. It’s definitely not efficient. Also the isotropic TV may not be the best choice because it seems to cause motion artefacts as it smooths the temporal evolution.
Nevertheless it is a fun experiment and illustrates the capabilities of CIL.
For another example of dynamic reconstructions see the demo on dynamic CT: https://github.com/TomographicImaging/CIL-Demos/blob/main/demos/3_Multichannel/02_Dynamic_CT.ipynb
[ ]: