About –– https://hackmd.io/@ccpi/cil-online24

The April 2024 CIL online training will take place on Monday 29th and Tuesday 30th April from 1-5pm UK time.

This document is for participants to find information about the course, ask questions and share answers with each other.

Resources

• CIL on GitHub incl. installation instructions: https://github.com/tomographicImaging/cil

• Documentation: https://tomographicimaging.github.io/CIL/v23.1.0/index.html

• CIL demos: https://github.com/TomographicImaging/CIL-Demos

• CIL User Show-case notebook collection: https://github.com/TomographicImaging/CIL-User-Showcase/

• Join Discord CIL support platform: https://discord.com/invite/9NTWu9MEGq

• Publications on CIL: https://doi.org/10.1098/rsta.2020.0192 and
  https://doi.org/10.1098/rsta.2020.0193

• Recent paper from the CIL team on a directional reconstruction method (using the PDHG optimization algorithm) for limited-angle CT, which won third prize in the Helsinki Tomography Challenge 2022: https://www.aimsciences.org/article/doi/10.3934/ammc.2023011

• Slides from today: Emailed to participants.

Training program

Monday 29th April

Time	Topic	Presenter
1:00 - 2:30	Geometry, data reading and FBP/FDK reconstruction Overview of the day Icebreaker Introduction to CT and FBP Introduction to CIL Notebook demonstration: 1_Introduction/01_intro_walnut_conebeam.ipynb STFC cloud setup Break out rooms for notebook: Go to: CIL-Demos/examples/1_Introduction/01_intro_walnut_conebeam.ipynb	Jakob Sauer Jørgensen and Franck Vidal
2:30 - 2:45	Break
2:45 - 3:45	Pre-processing Introduction to pre-processing Break out rooms for notebook: Go to: CIL-Demos/demos/1_Introduction/02_intro_sandstone_parallel_roi.ipynb	Jakob Sauer Jørgensen and Franck Vidal
3:45 - 4:00	Break
4:00 - 4:50	Time to explore and discuss Introduction Option 1 in breakout rooms: CIL-Demos/demos/1_introduction/03_preprocessing.ipynb Option 2 in breakout rooms: CIL-Demos/demos/1_introduction/exercises/02_preprocessing_seeds_conebeam.ipynb Option 3 in breakout rooms: CIL-Demos/demos/1_introduction/exercises/03_where_is_my_reader.ipynb	Jakob Sauer Jørgensen and Franck Vidal
4:50 - 5:00	Wrap up and where to go from here	Jakob Sauer Jørgensen and Franck Vidal

Tuesday 30th April

Time	Topic	Presenter
1:00 - 2:15	Introduction to iterative reconstruction Overview of the day Icebreaker Introduction to iterative reconstruction Breakout rooms for notebook: CIL-Demos/demos/1_Introduction/04_FBP_CGLS_SIRT.ipynb	Edoardo Pasca and Margaret Duff
2:15 - 2:30	Break
2:30 - 3:45	Regularised Iterative Reconstruction Introduction to regularisation Breakout rooms for notebook : CIL-Demos/demos/2_Iterative/01_optimisation_gd_fista.ipynb	Edoardo Pasca and Margaret Duff
3:45 - 4:00	Break
4:00 - 4:45	Time to explore and discuss Introduction Option 1 in breakout rooms: CIL-Demos/demos/1_Introduction/05_usb_limited_angle_fbp_sirt.ipynb Option 2 in breakout rooms: CIL-Demos/demos/2_Iterative/05_Laminography_with_TV.ipynb Option 3 in breakout rooms: CIL-Demos/demos/3_Multichannel/03_Hyperspectral_reconstruction.ipynb	Edoardo Pasca and Margaret Duff
4:45 - 5:00	Next steps and questions	Edoardo Pasca and Margaret Duff

Exercises

Icebreaker exercise - Monday

In your breakout rooms, please:
-Introduce yourselves and your
background
-Explain what you want to get out of the training
-Tell the group about the strangest dataset you have ever come across
-Trainers will be in main room for support

1. What is the strangest dataset you have ever come across? Please write your answer below:

• Pokeballs!
• Kinder egg
• Grand theft auto 5 for pose tracking
• Medieval swords
• A model train I built
• Worm burrows
• Fruit fly
• Fairy trapped in a slime jar
• Tennis ball
• teeth of snail
• Daphnia
• My friend's skull (he's alive)
• Solidification metal dendrites
• Ancient cylindrical seal
• fossil human ancestor
• Lego man!
• Mice skull (with a hole on it)
• Drug Pill
• Lego gandalf
• dens invaginatus

2. What are you hoping to get out of this course on CIL? Write your answer below:

• New reconstruction methods
• Conjugate gradient methods in reconstruction
• More informed knowledge of reconstruction & what software options are available for reconstruction
• Flexibility in reconstruction methods with limited data++
• Understand the basics of tomo reconstruction with CIL
• Better understanding of iterative reconstruction methods using CIL
• Better understanding of CIL to reconstruct with very limited datasets
• Understanding how to use for larger then memory datasets
• Getting better knowledge about CT and reconstructions using CIL
• new reconstruction methods
• spectral/hyperspectral imaging and acquisition

[]
[]

Icebreaker questions - Tuesday

Can you guess what these objects are

Don't write your answers here, give other people a chance to guess!

Questions

What questions do you have?

Write your questions here and we will answer below

Introduction to iterative

What is the objective function?

The objective function is the thing that you want to minimise it usually contains a data fitting term which measures how well your solution matches the measured data. It sometimes also contains a regularisation term which can provide extra information

How to check if the minimisation fuction is in a local minimum or it has a monotonic trend? or How can we know the algorithm has reached the global minimum rather than a local minimim?

That's a good question. If the operator is linear (like a tomographic operator) and the objective function is convex then if we are at a local minima then it must be global. We assume we are at a minima if our optimisation algorithm has converged, if in each iteration very little changes.

In the case of a non-linear operator or non-convex objective function, local minima are possible and optimisation is difficult. We could try different initialisations or algorithms with momentum to avoid local minima

Can you please explain why the step-size=None in the initialization of the gradient descent?

Excellent question. CIL's gradient descent automatically does a backtracking line search to find a suitable step-size if one is not provided. This guarantees convergence. Instead, if you want to pass a step-size, convergence is not always guaranteed e.g. if you pick a step-size too large you could keep jumping over the minimum point.

A step size of 1/(2*Lipschitz constant of the objective function) is usually a good number if you want a constant step size.

In example 4 I reduced the number of projections to 10, but then gave both algorithms the actual phantom as the initial guess. SIRT was able to basically predict the phantom, whereas CGLS still shows the reduced projection artificacts. Why is that?

An interesting question - with the noise in the image, the |Au-b|_2^2 is not zero even when u is the ground truth thus the algorithms may move away from the initial ground truth point. CGLS generally converges faster than SIRT so the results from SIRT may look better because it has not yet converged. I wonder also if you are running SIRT with a non-negativity constraint which would also improve your SIRT solution.

Can we use custom stop criteria for iterative recon?

Yes! Take a look at our new callbacks https://tomographicimaging.github.io/CIL/nightly/optimisation/?highlight=callbacks#utilities where you can pass a custom stop criteria

Can iterative methods be used to mitigate artefacts due to improper rotation and ring artefacts from synchrotron sources

We suggest to use the appropriate methods to correct for this type of artifact: centre of rotation processor and ring-remover in CIL or methods found in other packages like tomopyor algotom.
After an appropriate pre-processing has been done, iterative methods may be used and further reduce the artifacts.

More iterative reconstruction and regularisation

Can we do iterative reconstruction in a batch mode, after finding the optimal recon parameters in the first dataset? And I assume the following reconstructions after the 1st will be faster? I mean, will the iteration process be faster after the 1st one? i.e., learn from the 1st data.

If there are similar characteristics between the datasets then you can probably use the same parameters. We don't have a mechanism to batch process in CIL at the moment but you can use CIL if you parallelise the code yourself.

One of the major advantages of the iterative reconstruction is to make the in-situ imaging even faster, with reduced number of projections and faster exposure. So, I am not sure if any one has tried to install CIL into any synchrotron beamline, which can help to optimise the scanning parameters from the beginning? Not sure if it is allowed by the beamline scientist though.

We are currently working with Diamond, DIAD and I12 beamlines and with the ID15 beamline at the ESRF. We are happy to discuss further and support you in this work.

Hi, how many doppelgangers can you populate this landscape with so they can take perhaps many paths down the mist simultaneously? Are doppelgangers a complete no go?

Swarm type of algorithms are useful when the objective function is not convex. For our XCT problem the objective functions are convex so we do not have swarm algorithms.

Can you please touch on the criterion to choose between the analytical and the stochastic algorithm for object reconstruction? Thank you

This depends on the application. As a rule of thumb, you should use analytical methods, and only if this leads to unsatisfactory results, you should investigate iterative methods.

Choosing the regularisation parameter: Is there an upper/lower limit for the search?

Unfortunately, it is hard to set limits a priori. For least squares data fitting term and TV regularisation, I observed that a ball park for

$\lambda$ is by taking the 10% of the ratio of the norm of the projection operator and the norm of the GradientOperator.

Can you please explain the contribution of PDHG to Total Variation?

PDHG is a useful algorithm for optimisation objectives where both the data discrepancy and the regularisation term are not differentiable. You can also build up quite complex objectives, such as seen in this recent CIL paper: https://www.aimsciences.org/article/doi/10.3934/ammc.2023011/. In the case of just TV reconstruction of a least squares problem, you could use either FISTA or PDHG as optimisation methods with similar results.

In notebook 01_optimization_gd_fista. I get an errorat step 47 even if alpha = 0.02 end number iterations = 300

Make sure you run the algorithm for the right number of iterations with the cell myFISTATV.run(200,verbose=2) and if alpha is set to 0.02 the test should pass

Your reconstruction should look like this:

If it looks right just move on to the next cell and don't worry about the check.

Laminography

In laminography do we have 2 dimensions of same order and the 3rd quite small? Does it also include the case of a rectangular parallelopiped where all 3 dimensions a quite different from one another?

A typical laminography cone beam configuration might look like this:

If your object has very different dimensions you might have something like this:

At low tilt angles you will still get a range of different pathlengths, however as you tilt more they will become closer to your z-dimension thickness. As you tilt more you lose information at a wider range of angles so it's worth simulating your sample and seeing what works with a balance between pathlengths (and contrast range) and missing data.

Hyperspectral reconstruction

Is the case that each energy channels is recontructed indepedently in the example? There might be some cases where the solution is one channel might contribute to the solution in a different channel.

Absolutely. As the notebook progresses, more complex regularisation will be added, ending up with TV on space and TGV on regularisation on the channel dimension.

Can we do iterative reconstruction in a batch mode on CIL, after finding the optimal recon parameters in the first dataset? And I assume the following reconstructions after the 1st will be faster?

Introduction to CT and FBP

I have my own data which is saved as a numpy array, how do I reconstruct it with algorithms like CGLS?

Thanks for your question! Later today you can have a look at the "Where's my reader" notebook which discusses loading your own data and sorting out your geometries. Then tomorrow we will look at iterative algorithms like CGLS.

"Novel imaging modalities" - does this include "Edge Illumination phase contrast imaging"?

Thanks, we are currently actively working at implementing a phase retrieval method - not currently for edge illumination methods but we would be interested to discuss!

Are the any phase retrieval methods (e.g. Paganin's retrieval) implemented in CIL?

Thanks, we are currently actively working at implementing a phase retrieval method.

When doing a FBP reconstruction, how can we decide which filter is best to use, ie. ram-lak, shepp-logan etc..?

There is generally a trade off between more blurring and more noise reduction. It is often trial and error for each dataset based on what you need the data for.
You can have a go with the walnut notebook, here in 2D for ease of computation.

data2d = data.get_slice(vertical='centre')
data2d.reorder('tigre')
ig2d = data2d.geometry.get_ImageGeometry()
fdk =  FDK(data2d, ig2d)
recon = {}
for filt in ['ram-lak', 'shepp-logan', 'cosine', 'hamming', 'hann']:
    fdk.set_filter(filt)
    recon[filt] = fdk.run()
show2D([r.array[400:600, 400:600] for r in recon.values()], title=[r for r in recon.keys()], fix_range=(-0.01, 0.06), num_cols=3, size=(8,5))

I produced the following plot:

In islicer, what does the minmax range refer to ?

islicer(recon, direction='vertical', minmax=(-0.01, 0.06)
minmax sets the initial range of the colorbar, you can check any other parameters in our documentation https://tomographicimaging.github.io/CIL/nightly/utilities/#islicer-interactive-display-of-2d-slices
Unfortunately, this is not consistent with show2D where the same parameter is called fix_range

Introduction to pre-processing

Write your questions here and we will answer below

DataReader - how to select the right one. We have a Phoenix Vitomex s160 in the lab. Do CT machine manufactureres adher to a standard for data organization.

Unfortunately each manufacturer chooses their own data format. We can work with you to create the appropriate reader for your machine. What will need to be done is: 1) to find the geometry information and 2) the projection data. Later on this afternoon you will have an option to look at the "Where's my reader" notebook which looks at setting up geometries for custom data. Similarly, the pre-processing sandstone notebook creates a geometry for a dataset stored as a mat file.

Are "dark" and "flat" data for pre-processing stored in the specific data structure used in the example?

Thanks, this will depend on the actual scanner. We mainly support NIKON and ZEISS: NIKON stores normalised data, ZEISS stores the flat fields and our reader will apply the normalisation. The other readers we have are NEXUS data reader, which is a reader for data stored in HDF5. Raw and TIFF readers do not have any information about flats and darks. Check here for a discussion on how to store this info

Source/Ilumination: Beyond Point source? Anyway to enforce finite source size? Any method to create structured illumination e.g. mask between source and sample.

We do not have any way to do this natively, although I can think of methods to achieve it. We've had examples of creating a mask between the source and sample in an algorithm by implementing a MaskOperator (https://tomographicimaging.github.io/CIL/nightly/optimisation/?highlight=maskoperator#cil.optimisation.operators.MaskOperator)

Does CIL have pre-processors for stitching together data e.g. regions of interest?

Not at the moment - do you have an examples of other software or other places where it is implemented? Algotom has a function for specific types of stitching (https://algotom.github.io/toc/section1/section1_4.html#sinogram-stitching-for-a-half-acquisition-scan) - is this the kind of thing you mean? We have a description of a similar method in our user showcase notebooks (https://github.com/TomographicImaging/CIL-User-Showcase/blob/main/009_offset_CT_apple/Apple_offset_reconstruction.ipynb)

Are there any CIL functions to correct for sample drifts during CT scans

Not out of the box in CIL but this is a very interesting question. We would be interested in talking more about this

Does CIL Work with non conventional rotation paterns?

Non standard trajectories that can be described with individual positions of source, panel and object are currently being worked on.

Is there any benefit to using astra vs tigre?

They have different licenses, ASTRA is GPL and TIGRE BSD so this might be one consideration for you. TIGRE have papers that references the projectors they have implemented, whereas ASTRA does not as far as I know – this may be something that’s important to you if you want to understand the projectors. Otherwise both should work with most set-ups/geometries. ASTRA may have a slightly smaller peak RAM use via the CIL wrappers, but this is something we are working to address with TIGRE.

I cannot understand how looking at the sinogram we can figure out what is wrong with the image reconstruction

Thanks for your question, we answered this very quickly in the slides. A lot of this is experience but there are some resources online, that you can google. We can also share the slides after the training

How many channels does cil support?

We can support multi-channel images e.g. hyperspectral or dynamic CT, with any number of channels. There will be an opportunity to see this at the end of the day tomorrow
CIL Data containers have a maximum of 4 dimensions, 3 spatial + 1 channel dimension which could be energy or time, for instance.

Can the acquired ROI be cropped instead of padding data?

If we crop the acquisition data we will be throwing away data and this is undesired

Can you please review the ROI problem? Is the FOV smaller than the object we want to image or is the object much smaller than the FOV? In the latter case I understand the padding. In the first case I don't.

We consider the former case, the FOV is smaller than the the whole object. We pad to correct for the missing data. . This is because the data is region-of-
interest, i.e., the sample was larger than the field of view, so projections are truncated (i.e. do not have air on both sides as also seen in the sinograms). A simple way to compensate for this is to extend or pad the data on both sides of the projections.

The following plot depicts the scan configuration for the walnut dataset. If we move the object closer to the source (increasing the magnification) but we keep the size of the reconstruction volume constant, the object is larger than the detector. The X-rays go through the material at the corners of the object only at certain rotation angles, creating problems to algorithms. The last plot shows that by padding the data effectively elongating the detector, the object remains in view.

Padding the projections will help if the sample outside the FOV is a low/constant density. Highly attenuating features outside the FOV will cause issues with the reconstruction as the padded data is not a good approximation.

How does the centre of rotation correction take sharpness into account?

With 360 degrees if the geometry doesn’t describe the data well the features in the reconstructed image look blurred. By changing the centre-of-rotation offset and measuring the sharpness of the image we can find the geometry that best describes the data

The image-sharpness algorithm backprojects the sobel-filtered projections with different offsets. This means we end up with strong edges in the case of aligned geometry and blurring otherwise. We maximise the autocorrelation of the reconstruction.

Can the centre of rotation correction find the roll/tilt of the rotation axis too?

The centre of rotation correctors only finds the offset in one slice (so a simple offset). However, you could run it twice and find the offset in two slices, then manually update the geometry using:

• my_acquisition_geometry.set_centre_of_rotation(offset, angle)

• my_acquisition_geometry.set_centre_of_rotation_by_slice( offset1, slice_index1, None, slice_index2)

The current centre of rotation corrector algorithms can only be applied to un-tilted geometry. This is something we hope to remove the restriction on in the future.
I notice this is missing from our documentation, so we’ll update it to be clearer.

About Hyperspectral recontruction, do the algorithms work for a source different from an microfocus X-ray tube? Example: a light source?

Yes, you can use CIL for light tomography and indeed Franck has a project with a PhD student in which they had a play with a light tomography set up. We have also used CIL with neutron tomography data.

About the optimization methods, there are methods like CGLS, GD and so on in CIL. are any of them suitable for a nonlinear problem? in case we define our regularization term in a way that leads to a nonlinear problem.

If you can define an adjoint for non-linear operator then the CIL algorithms will work. However, with the non-linearity you loose the guarantanees for convergence. Basically, you can give it a go but all bets are off!

Got a deprecation warning.

Thank-you for letting us know!

Jupyterhub Instructions

• Go to:  https://tinyurl.com/cil24online write your name next to a username to claim it for the exercises
• CIL Jupyter notebook server: https://training.jupyter.stfc.ac.uk/
• Sign up with the username you claimed and a password of your choice. Create a memorable password. We can't change it.

FAQ Jupyterhub

I cannot create a log in into jupyterhub

Make sure to have chosen a user name from the online form circulated at the beginning of the training.

Did you use a username of the form cil24-online-XX?

The authorised user names on the platform are cil24-online-xx where xx is a number between 1 and 75, that is cil24-online-1, …, cil24-online-75. If you used any different user name you won't get access to jupyterhub.

Did you forget the password?

We do not have any way to recover the user password, so please remember it/write it down. Should you forget, please get in touch with one of the trainers.

403 : Forbidden. XSRF cookie does not match POST argument . This is the issue I am getting

Please try and refresh your page, if you are still getting this error please get back in touch

Feedback

Feedback Tuesday

Please fill in this quick microsoft form:
https://forms.office.com/e/WtKHKyKhZY

Monday

Please share one thing you liked about the course today

• Liked the step by step breakup . Easy to understand

Liked the step by step breakup. Easy to Understand

• Awesome explanation on every steps
• The notebooks are very well commented and easy to follow
• The interactive experience of the notebooks
• Useful and easy to understand examples in the notebooks
• easy to understand , need it to be a little longer
• The explanations are clear as well as the mathematical models and the codes, thank you!
• Thank you! great workshop.
• Very interesting real world datasets and problems
• Exercises were easy to follow

Please share one thing that could be improved about the course today

• Deeper overview of the documentation - processors vs frameworks vs recon
• Maybe some bibliography, but there is lot of material

Maybe some bibliography, but there is lot of material