<jupyter_start><jupyter_text>Using Scipy Solvers for MR Reconstruction with SIRFThis demo is a 'script', i.e. intended to be run step by step in aPython notebook such as Jupyter. It is organised in 'cells'. Jupyter displays thesecells nicely and allows you to run each cell on its own. Learning objectiveThis notebook is intended to give you the basis to implement iterative MR reconstruction algorithms using SIRF.By the end of this notebook you will know how to- define an objective function and its gradient using the SIRF MR acquisition models- call a scipy-solver to optimise the objective function.Of course, you can use CIL optimisers as opposed to scipy, but that is for another notebook. First version: 21st of May 2021Author: <NAME>CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  Copyright 2015 - 2021 Rutherford Appleton Laboratory STFC.  Copyright 2015 - 2021 University College London.  Copyright 2015 - 2021 Physikalisch-Technische Bundesanstalt.This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (http://www.ccpsynerbi.ac.uk/).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.0Unless required by applicable law or agreed to in writing, softwaredistributed 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 andlimitations under the License.<jupyter_code>#%% make sure figures appears inline and animations works
%matplotlib notebook

# Setup the working directory for the notebook
import notebook_setup
from sirf_exercises import cd_to_working_dir
cd_to_working_dir('MR', 'e_advanced_recon')
__version__ = '0.1.1'

# import engine module
import sirf.Gadgetron as pMR

# import further modules
import os, numpy as np

import matplotlib.pyplot as plt
import matplotlib.animation as animation

from IPython.display import HTML
from sirf_exercises import exercises_data_path
data_path = exercises_data_path('MR', 'PTB_ACRPhantom_GRAPPA')<jupyter_output><empty_output><jupyter_text>Auxiliary functionsTo facilitate plotting<jupyter_code># This is just an auxiliary function
def norm_array( arr ):
    arr = np.squeeze(arr)
    min_a = abs(arr).min()
    max_a = abs(arr).max()
    if (max_a-min_a) < np.finfo(np.float32).eps:
        return arr
    else:
        return (arr - min_a)/(max_a - min_a)

# LOADING AND PREPROCESSING DATA FOR THIS SET
filename_grappa_file = os.path.join(data_path, 'ptb_resolutionphantom_GRAPPA4_ismrmrd.h5')
acq_data = pMR.AcquisitionData(filename_grappa_file)
preprocessed_data = pMR.preprocess_acquisition_data(acq_data)
preprocessed_data.sort()
print("Number of acquisitions is {}".format(preprocessed_data.number()))
# WE DO A GRAPPA RECONSTRUCTION USING SIRF

recon = pMR.CartesianGRAPPAReconstructor()
recon.set_input(preprocessed_data)
recon.compute_gfactors(False)
print('---\n reconstructing...')

recon.process()
# for undersampled acquisition data GRAPPA computes Gfactor images
# in addition to reconstructed ones
grappa_images = recon.get_output()<jupyter_output><empty_output><jupyter_text>Generation of an aquisition modelAs before we want to construct an `AcquisitionModel` which has methods for both the simulation of the acquisition process, as well as its adjoint operation.<jupyter_code># First we generate different channels
csm = pMR.CoilSensitivityData()
csm.smoothness = 10
csm.calculate(preprocessed_data)
# Now we have enough information to set up our model
E = pMR.AcquisitionModel(preprocessed_data, grappa_images)
E.set_coil_sensitivity_maps(csm)

# Now we can hop back from k-space into image space in just one line:
bwd_img = E.backward( preprocessed_data )

# PLOT THE RESULTS
grappa_img_arr = norm_array(grappa_images.as_array())
bwd_img_arr = norm_array(bwd_img.as_array())

fig = plt.figure(figsize=(9, 4))
plt.set_cmap('gray')

ax = fig.add_subplot(1,2,1)
ax.imshow(abs(grappa_img_arr), vmin=0, vmax=1)
ax.set_title('Result of GRAPPA Reconstruction')
ax.axis('off')

ax = fig.add_subplot(1,2,2)
ax.imshow(abs(bwd_img_arr), vmin=0, vmax=1)
ax.set_title('Result of AcquisitionModel.backward()')
ax.axis('off')

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>With the AcquisitionModel we can use solvers provided by PythonWe want to use existing implementations of optimisation algorithms in `scipy.optimize.minimize` to solve our reconstruction problem:<jupyter_code>%%html
<iframe src="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html" width="1000" height="800"></iframe><jupyter_output><empty_output><jupyter_text>Take-home messages from reading their documentation:- `scipy.optimize.minimize` will minimize the function `fun(x, *args)` with respect to  `x`, i.e. find a solution to the problem $\hat{x} = \min_{x} C(x)$- `x = array_like, shape (n,)` is REAL valued (which is a bit of a bummer because our beautiful MR images are complex valued)We need to give it a list of arguments- You need to supply the objective function `fun(x,*args) -> float`.- You need to at least supply the gradient of the objective function `jac(x, *args) -> array_like, shape (n,)`- Since `x` is real data we need to split our complex images into two channels: `[img] = [real(img), imag(img)]`Our cost function is:- $C(x) = \frac{1}{2}\lVert E x - y \rVert_2$ - $dC(x) = E^H (Ex - y)$    TASK: Program the splitting of real and imaginary partSetup: The AcquisitionModel provides SIRF objects. We get and set the data using the methods `as_array()` and `fill()`.Problem: We need to split the data into real and imaginary part, s.t. we only supply real-valued arrays to `scipy`.Task: write the code for the two functions in the cell below to split complex valued arrays into one twice-as-longHint 1: this can be achieved using the methods `np.concatenate`, `np.split`, `np.real`, `np.imag`, and the `array.flatten()` methods.Hint 2: All arrays need to be flattened, i.e. of shape `(n,)` since that's what `scipy` requires. So no need to worry about shapes, when we use `fill()` SIRF will track the shape for us.<jupyter_code># WRITE YOUR FUNCTIONS HERE:
import scipy.optimize as optimize

def split_complex_to_real(data):
    pass # replace the word pass with your code


def merge_real_to_complex(data):
    pass # replace the word pass with your code

# SOLUTION CELL, TRY IT FIRST ON YOUR OWN THOUGH


































import scipy.optimize as optimize

def split_complex_to_real(data):
    
    data = np.concatenate( (np.real(data.flatten()), np.imag(data.flatten())))
    return data.astype(np.float64)

def merge_real_to_complex(data):
    data = data.astype(np.float64)
    data = np.split(data,2)
    return data[0] + 1j*data[1]<jupyter_output><empty_output><jupyter_text>TASK: Program the splitting of real and imaginary part Setup: The objective function `fun(x, *args)` takes `x` as the argument, and all the other arguments are passed in `*args` Problem: None so far. Task:Write the code for the cost function and it's graident in the cells below. Maybe, first, think of all variables except for `x` that are involved in the computation of `C(x)` and `dC(x)` and pass them using `args`. Hints:- Check the minimal example below to familiarise yourself with the way how to pass `args` to a function. The objects in the `args` tuple can be of any type.- keep in mind, that `scipy` will use  the same variable for `args` in both cost function and gradient.- make sure to use the splitting and merging of complex data inside the functions.- The acquisition model needs SIRF objects, so you need to use `fill()` and `as_array()` to get and set the data.<jupyter_code># Minimal example using the *args syntax
num_fruit = 15
name_fruit = "Chicken"
x = np.zeros((3,3))

def fun(x, *args):
    
    fruit_number = args[0]
    fruit_type = args[1]
    
    message = \
    "Our x is of type {}, and shape {}. We have {} things that are called <{}>.\n " \
    "The objects don't need to be strings and integers, but can be anything, e.g. an acquisition model or rawdata.'"\
    .format(type(x), x.shape, fruit_number, fruit_type) 
    
    print(message)                    

extra_arguments = (num_fruit, name_fruit)        
fun(x, *extra_arguments)
        
def objective(x, *args):
    assert len(args)==3, "Please give three arguments only"
    pass # replace pass by your code

def grad_objective(x, *args):
    
    assert len(args)==3, "Please give three arguments only"
    pass # replace pass by your code
# SOLUTION CELL, TRY IT FIRST ON YOUR OWN THOUGH











































def objective(x, *args):
    
    
    assert len(args)==3, "Please give three arguments only"
    E = args[0]
    i = args[1]
    y = args[2]
    
    data = merge_real_to_complex(x)
    i = i.fill(data)
    
    c =  E.forward(i)-y
    return 0.5 * c.norm() ** 2

def grad_objective(x, *args):
    
    assert len(args)==3, "Please give three arguments only"
    E = args[0]
    i = args[1]
    y = args[2]
    
    data = merge_real_to_complex(x)
    i = i.fill(data)    
    
    dObj = E.backward(E.forward(i) - y)
    dObj = dObj.as_array().flatten()
    
    return split_complex_to_real(dObj)
# fix some extra arguments

# normalising the rawdata
y = preprocessed_data
y = np.sqrt(2) * y / y.norm() # this means that C(x=0) = 1

print("The norm of our k-space data is: {}".format(y.norm()))

template_img = grappa_images.copy()

extra_args = (E, template_img, y)<jupyter_output><empty_output><jupyter_text>TASK: Check if your gradient is correctProblem: We don't know if the objective function and its gradient are correct.Setup: But we do know that:- we normed our rawdata s.t. $C(x=0)=1$ - we can see that such that $dC(x=0) =E^H (E(0) - y) = -E^H y$. That means, the gradient of the objective at zero-image is the backward of the rawdata, i.e. the zero-filled reconstruction)Task: evaluate both Hint: keep in mind that both functions expect a real-valued array.<jupyter_code># set the start value
img = grappa_images.copy()
img.fill(0.0 + 1j*0.0);
# call the correct methods
x0 = pass # keep it real!
grad0 = pass 
x0 = pass
#SOLUTION CELL: DON'T RUN IF YOU DIDN'T TRY!











































# compute the gradient at img=0 and initialize the reconstruction with that image
x0 = split_complex_to_real(img.as_array().flatten().copy())
grad0 = grad_objective(x0, *extra_args)
x0 = grad0 
c0 = objective(x0, *extra_args)
print("Objective function start is: {}".format((c0)))

# Plotting the result
img.fill(-1*x0) # the minus sign is important, hopefully ever

fig = plt.figure(figsize=(9, 4))
plt.set_cmap('gray')

ax = fig.add_subplot(1,1,1)
ax.imshow(norm_array(np.abs(img.as_array())), vmin=0, vmax=1)
ax.set_title('Gradient of the cost function at x0=0')
ax.axis('off')
# run the scipy optimiser
cg_tol = 0 
cg_iter = 10
extra_options = {"maxiter": cg_iter, "disp": True, "gtol":cg_tol} # just some options to give to the optimiser

res = optimize.minimize(objective, x0, args=extra_args, method='CG', jac=grad_objective, options=extra_options)
img = img.fill(res.x)<jupyter_output><empty_output><jupyter_text>Let's display the result<jupyter_code>fig = plt.figure(figsize=(9, 4))
plt.set_cmap('gray')

ax = fig.add_subplot(1,1,1)
ax.imshow(norm_array(np.abs(img.as_array())), vmin=0, vmax=1)
ax.set_title('Result of the optimisation')
ax.axis('off')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1002751977984176, 'fraction_numerical': 0.012039903680770554, 'mean_word_length': 4.421445221445222, 'pattern_counts': {'":': 3, '<': 38, '<?xml version=': 0, '>': 39, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '220807', 'license': 'Apache-2.0', 'n_tokens_mistral': 3869, 'n_tokens_neox': 3696, 'n_words': 1322, 'path': 'notebooks/MR/e_advanced_recon.ipynb', 'repo_name': 'DANAJK/SIRF-Exercises', 'size': '11636'}