ZeissDataReader Shifts and Metadata

epapoutsellis commented 12 months ago

For about one month, I have been trying to reconstruct a dataset (cone-beam geometry) using the ZeissDataReader, and I found some bugs in terms of reading the actual data and reading the metadata.

First of all, using the dxchange instead of xrmreader, we lose some important information for the metadata. For my data, if I use dxchange.read_txrm(filename_sinogram) function the keys of the metadata dictionary are

dict_keys(['facility', 'image_width', 'image_height', 'data_type', 'number_of_images', 'pixel_size', 'reference_filename', 
'reference_data_type', 'thetas', 'x_positions', 'y_positions', 'z_positions', 'x-shifts', 'y-shifts', 'reference'])

If I use the xrmreader.read_metadata(filename_sinogram) function the keys of the metadata dictionary are

dict_keys(['facility', 'image_width', 'image_height', 'data_type', 'number_of_images', 'pixel_size', 'reference_exposure_time', 
'reference_current', 'reference_voltage', 'reference_data_type', 'image_data_type', 'thetas', 'x-shifts', 'y-shifts', 'AMC-x-shifts', 
'AMC-y-shifts', 'temperature-x-shifts', 'temperature-y-shifts', 'source-x-shifts', 'source-y-shifts', 'align-mode', 'center_shift', 
'rotation_angle', 'source_isocenter_distance', 'detector_isocenter_distance', 'cone_angle', 'fan_angle', 'camera_offset', 'source_drift', 
'current', 'voltage', 'power', 'exposure_time', 'binning', 'filter', 'scaling_min', 'scaling_max', 'reference’])

I do not know why, but the latter seems more complete. Besides the x-shifts, y-shifts and the reference we have more important shifts such as AMC, temperature and source. For more information, see here. There are also some other keys which I have no idea what they describe, such as align-mode, cone_angle, fan_angle, camera_offset, scaling_min, scaling_max. If anyone is familiar with these keys, I would really appreciate some explanation.

Example

If we use the ZeissDataReader to read the WalnutData, everything is ok. Using the xrmreader.read_metadata(filename_sinogram) function, the AMC, temperature are None but the source-x-shifts and source-y-shifts are close to 0, meaning that we do not need to include this correction.

AMC-x-shifts =  None
AMC-y-shifts =  None
temperature-x-shifts =  None
temperature-y-shifts =  None
source-x-shifts =  [0.00010398 0.00020615 0.00030651 0.00040505 0.00050179 0.0005967
 0.00068981 0.0007811  0.00087058 0.00095825]
source-y-shifts =  [-1.95432222e-05 -3.88436601e-05 -5.79013140e-05 -7.67161837e-05
 -9.52884438e-05 -1.13617745e-04 -1.31704321e-04 -1.49548112e-04
 -1.67149119e-04 -1.84507342e-04]
align-mode =  None
center_shift =  None
rotation_angle =  None

However, in my data I have something like:

x-shifts =  [-0.02338324 -8.31449032  3.44475412  7.60942364  4.46493196  1.39798403
  6.30554199 -0.50505334 -2.10339952  9.3211956 ]
y-shifts =  [-0.06622199 -0.95438933  2.88818741  0.95634484 -2.7804389   0.52345252
 -1.7584523  -1.91345716  2.86733198 -0.94697839]
AMC-x-shifts =  [-0.03012854 -0.02821015 -0.02629176 -0.02437337 -0.02245498 -0.02053659
 -0.0186182  -0.01669981 -0.01478142 -0.01286303]
AMC-y-shifts =  [0.10338005 0.10249254 0.10160504 0.10071754 0.09983003 0.09894253
 0.09805502 0.09716751 0.09628001 0.09539251]
temperature-x-shifts =  [-0.07080591 -0.06542346 -0.07087189 -0.0660259  -0.07146409 -0.07147949
 -0.06607211 -0.07149497 -0.07149501 -0.06607226]
temperature-y-shifts =  [-0.11793955 -0.11193907 -0.11793955 -0.10978469 -0.11578517 -0.11578517
 -0.10978469 -0.11578517 -0.11578517 -0.10978469]
source-x-shifts =  [0.00125837 0.00251942 0.00378316 0.00504959 0.00631873 0.00759054
 0.00886504 0.01014224 0.01142212 0.0127047 ]
source-y-shifts =  [0.00138348 0.00277969 0.00418863 0.00561029 0.00704469 0.00849178
 0.00995162 0.01142418 0.01290946 0.01440747]
align-mode =  [3]
center_shift =  0.8091099262237549
rotation_angle =  0.0

meaning that somehow we need to correct our projections for all of these shifts.

Solution

I have tried a lot of things which I will not mention here, in order to understand what is the problem with the CIL reconstruction and why is so different from the recon.txrm which I guess is from a Zeiss software. The problem is in the following lines

https://github.com/TomographicImaging/CIL/blob/be57697c96d3e66032bac5de5246d0d276e3ba1e/Wrappers/Python/cil/io/ZEISSDataReader.py#L273-L276

We use np.roll which works only with integer values. Meaning that with an int rounding we do not perform a good interpolation. Actually, all the values in the x-shifts and y-shifts vectors that are in $\in(-1. , 1.)$ will end up to be 0 after the rounding. Below are two histogram plots of the y-shifts restricted to $\in(-1. , 1.)$ for the Walnut dataset with/without the int rounding

Why the Walnut recon looks correct using the np.roll and the `int` rounding???

Because the number of projections (108) where the y-shifts is $\in(-1. , 1.)$ is small compared to the total number of projections(800).

In my dataset the histogram plots are:

which means that the y-shifts for 300 projections out of 800 projections were approximated by 0.

Correct Implementation

is using the shift by scipy. This is how xrmreader implements the interpolation for all the shifts. Although, I am not sure if nearest mode and order 1 are good defaults.

def revert_shifts(projections, x_shifts, y_shifts):
    '''

    :param projections: projection data
    :param x_shifts: x shift parameter per projection
    :param y_shifts: y shift parameter per projection
    :return: shift corrected porjections
    '''
    num_projections, img_size_x, img_size_y = projections.shape
    shifts = np.zeros((num_projections, 2))
    shifts[:, 0] = x_shifts
    shifts[:, 1] = y_shifts

    # shift each projection
    for i in range(num_projections):
        # input is extended by repeating its last value to avoid zeros after cropping
        projections[i, :, :] = shift(projections[i, :, :], (shifts[i, 1], shifts[i, 0]), order=1, mode='nearest')

    return projections

In my case, I need also to correct for AMC, temperature, source shifts as well as to correct for the offset due to the center_shift. Below is a comparison, with all these different options a) using np.roll for x/y shifts, b) using scipy.shift for x/y shifts, c) using scipy.shift for all shifts & correcting the offset and d) the txrm.recon.

zeiss_compare

@antonyvam for future reference.

paskino commented 12 months ago

@epapoutsellis thanks for this detailed report.

It looks like there are 2 issues we should deal with:

extra fields that are not read in by the reader. I had a look at the xrmreader package and it seems it has 1 version only and the license is similar to MIT.
interpolating data for subpixel shifts. Our preferred path would be to handle this with a per-projection geometry setup, which is not currently available in CIL, but in development.

Maybe you can join the next developer's meeting or the user support.

epapoutsellis commented 12 months ago

@paskino Apologies but have meetings Tuesday and Friday. Will try to join on Friday. This is a recent package from Andreas Maier group.

What do you mean by per-projection geometry setup? I strongly suggest to use a parallel processing library to do all these shift corrections for all projections. Otherwise, it takes A LOT OF TIME. You know my favourite library for it.

paskino commented 12 months ago

I mean that we don't have to apply the shifts at all, just record them in the geometry of each projection. Then we pass this info to the engine which will do the projection.

Feel free to suggest another meeting time if those don't work for you.

paskino commented 11 months ago

@epapoutsellis can you list the additional parameters you require? These could be added to dxchange with a PR.

Thanks for sharing a ZEISS reference.

paskino commented 11 months ago

After discussion with @epapoutsellis @gfardell @effepivi we propose to tackle the reader issues in different ways:

remove np.roll in favour of scipy.shift with default to do same as what we do now (no backward incompatibility)
Proceed adding new metadata to dxchange, see above
The future reader will be able to return a per-projection geometry OR a shifted dataset with one single geometry

TomographicImaging / CIL