remove the artifacts introduced by changing flat fields, we need to make the correction dynamic as well.

$2$ Conventional flat field correction

Conventional flat field correction $($FCC$)$ works by exploiting the Lambert$-$Beer law, which gives a projection

of the X$-$ray signal attenuated by the material of the object we're scanning.

$p=I_0$exp$(-{\displaystyle \underset{\phantom{}}{\overset{\phantom{}}{}}}(l)dl)$

to get the attenuation from equation $1,$ we simply divide by the intensity $I_0$ and take the logarithm to get

rid of the exponential. To compensate for the dark current of the detector, we additionally subtract the dark

field image $d.$

${\displaystyle \underset{\phantom{}}{\overset{\phantom{}}{}}}(l)dl=-ln(\frac{p-d}{I_0-d})$

In practice $I_0$ and $d$ are mean images created from several scans to average out temporal differences. But

when the source and detector have instable behaviours over time, this can not adequately remove the arti$-$

facts.

$3$ FFC for dynamic processes

To correct for the changing conditions, we have to assign a changing $I_0$ to every projection. The corrected

or normalized projection for each angle is then given by

$n_j=\frac{p_j-(\frac{?}{\text{b}\text{a}\text{r}}(d))}{I_j-(\frac{?}{\text{b}\text{a}\text{r}}(d))},$

where $j$ is the index of the projection, $\frac{?}{\text{b}\text{a}\text{r}}(d)$ is the average dark field image and $I_j$ is the flat field corresponding

to the projection with the same index. Now we can certainly remove the object after every angle and take

a flat field, but that would severely influence the scanning quality. Ideally we want a flat field for every

projection, but we do not want to take intermediate flat fields. So we first make several scans to estimate

the fluctuations in the flat fields and then make our scan. We then try to estimate a flat field for each of the

projections out of the flat fields we took in the beginning. $4$ Eigen Flat Fields

Vincent Van Nieuwenhove developed a technique for that, which he named "Eigen Flatfield Correction"

$($EFF correction$)[$I$.$ The algorithm is as follows:

Compute eigen flatfield

$($a$)$ Prepare the flatfield matrix $F$ with size $NxM$ with $N$ the number of pixels in a single flatfield and

$M$ the number of flatfields. Then subtract the mean of the dark field images.

$($b$)$ Compute centered flat$-$field matrix $A$ by simply subtracting mean flatfield $\frac{?}{\text{b}\text{a}\text{r}}(f)$ from each individual

flatfield.

$($c$)$ Compute eigen vectors $v_i$ and eigenvalues ${}_i$ of the covariance matrix $C=A{A}^T.$ Refer to the

original paper $[$I$]$ to make this computation more efficient!

$($d$)$ Calculate standard deviation of the flatfield matrix A across several flatfields. Then draw $k$ sam$-$

ple matrices $S$ from a multivariate normal distribution with identical size as matrix A and that

standard deviation.

$($e$)$ Calculate eigenvalues of covariance of each matrix $S$ and compute the average over the $k$ repi$-$

tions

$($f$)$ Select the eigen flatfields that have the eigenvalues larger than $95\%$ of the computed eigenvalues

of the random matrices. You can do that by looking at the average value that each eigenvalue gets

for all the random matrices; if an eigen$-$value is above this average value plus $2$ times the stan$-$

dard deviation, then it an acceptable eigenvalue. At the end of this step, you will get $1$ principle

EFFs, which depends on the number of flatfields that meet the above criterion.

Estimate dynamic flatfields

$($a$)$ A single dynamic flatfield is generated by linearly combining above $k$ eigen flatfields with certain

weights. The optimization function is constructed using formula $(8)$ in the paper. It minimizes

the gradient magnitude of the flat fields.

$($b$)$ This optimization can be solved by using a customized cost function with the optimize $.$minimize

of the scipy Python package. For more information, please refer to the documentation.

$5$ Tasks

Write a Python script that uses the dataset on provided on Blackboard and imple$-$ ments the Eigen

Flat Field Correction described above. The dataset is composed by $($in this order$)20$ dark field images,

$300$ flat fields, $251$ projections, and $300$ flat fields. You can neglect the projection and some of the flat

field images if the memory requirement is too high for your laptop.

Tip: Make yourself a favour and define functions to implement a specific task$/$functionality$.$

Think of an approach to comp