hkaggarwal / modl

MoDL: Model-Based Deep Learning Architecture for Inverse Problems

MoDL

MoDL: Model Based Deep Learning Architecture for Inverse Problems

Reference paper:

MoDL: Model Based Deep Learning Architecture for Inverse Problems by H.K. Aggarwal, M.P Mani, and Mathews Jacob in IEEE Transactions on Medical Imaging, 2018

Link: https://arxiv.org/abs/1712.02862

IEEE Xplore: https://ieeexplore.ieee.org/document/8434321/

Presentation: https://github.com/hkaggarwal/modl/blob/master/MoDL_ppt.pdf

What this code do:

In the above paper, we propose a technique to combine the power of deep-learning with the model-based approaches. This code suggest how we can use a deep convolutional neural netwrok (CNN) as a regularizer to solve an optimization problem.

This code solves the following optimization problem:

 argmin_x ||Ax-b||_2^2 + ||x-Dw(x)||^2_2 

A can be any measurement operator. Here we consider parallel imaging problem in MRI where the A operator consists of undersampling mask, FFT, and coil sensitivity maps.

Dw(x): it represents the denoiser using a residual learning CNN.

Recursive MoDL architecture:

Main benefits of the MoDL:

1. One of the first deep model that works with parallel MRI data.
2. Can account for more general image forward models by using conjugate graident
3. Needs less training data because of weight sharing across MoDL iterations.

Output on the Brain data:

Output on the Knee data:

The output GIF is showing the reconstructions from 10 iterations of alternating minimization as described in the MoDL paper.

Dependencies

We have tested the code in Anaconda python 2.7 and 3.6. The code should work with Tensorflow-1.7 onwards. The dataset is in the hdf5 format. You may require to install hdf5 library in python. In Anaconda you can give following command conda install h5py

The training code requires tqdm library. It is a nice library that is helpful in tracking the training progress. It can be installed using: conda install tqdm

In addition, matplotlib is required to visualize the output images.

Brain Dataset

This git repository also includes a single image in the file demoImage.hdf5. The testing script tstDemo.py will use this image by default and does not require full data download for the testing purpose.

We have released the parallel imaging dataset used in this paper. You can download the full dataset from the below link:

Download Link : https://zenodo.org/records/6481291

You will need the file dataset.hdf5 to run the training code trn.py. You can download the dataset from the link provided above. Please ignore the future warning by python. You do not need to download the dataset.hdf5 for testing purpose.

This dataset consist of parallel magnetic resonance imaging (MRI) brain data of five human subjects. Four of which are used during training of the model and fifth subject is used during testing. Above link contain fully sampled preprocessed data in numpy format for both training and testing. We also provide coil-sensitity-maps (CSM) pre-computed using E-SPIRIT algorithm. Total file size is 3 GB and contains following arrays:

trnOrg: This is complex arrary of 256x232x360 containing 90 slices from each of the 4 training subjects. Each slice is of spatial dimension 256x232. This is the original fully sampled data.

trnCSM: This is a complex array of 256x232x12x360 representing coil sensitity maps (csm). Here 12 represent number of coils.

trnMask: This is the random undersampling mask to do 6-fold acceleration. We use different mask for different slices.

tstOrg,tstCSM, tstMask: These are similar arrays for testing purpose. There are total 164 testing images.

The undersampling mask, for both training and testing cases, is shared corresponding to 6-fold acceleration case.

How to run the code

First, ensure that Tensorflow 1.7 or higher version is installed and working with GPU. Second, just clone or download this reporsitory. The tstDemo.py file should run without any changes in the code. On the command prompt CD to this modl directory i.e. the directory containig tstDemo.py. Then you can run the test code using the command:

$python tstDemo.py from the command prompt.

Knee Dataset

The MoDL architecture can be trained using Knee datasets freely available from other sources. The directory knee_trained_MoDL contains a trained model on Knee data. Here we performed the training with a structured mask that is kept same for all the slices during training. The file knee_demo_data.h5 contains one raw image, mask, as we as coil sensitivity maps corresponding to a single slice from a particular subject for demo purpose only. Just use the command $python knee_demo_code.py to see the performance of the MoDL on Knee dataset.

Files description

The folder savedModels contain the learned tensorflow model parameters. tstDemo.py will use it to read the model and run on the demo image in the file demoImage.hdf5.

supportingFunctions.py: This file contain some supporting functions to calculate the time, PSNR, and read the dataset.

model.py: This file contain the code for creating the residual learning CNN model as well as the algorithm for conjugate-gradient on complex data.

trn.py: This is the training code

tstDemo.py: This is the testing code

Contact

The code is provided to support reproducible research. If the code is giving syntax error in your particular python configuration or some files are missing then you may open an issue or directly email me at [email protected]

modl's Issues

Will MoDL-MUSSELS open source?

I try to reduplicate MoDL-MUSSELS, but the network can't converge. Could you open the source?

Not able to download knee data

Dear Dr. Hemant,

I am not able to download the knee training data from the link given in the README.md file. Can you please rectify it ?

Sincerely
Aditya Rastogi

ValueError while runing trn.py

Hello sir, I am getting the following error from lines 147-151 of the trn.py file:

Successfully undersampled data
Traceback (most recent call last):
File "run.py", line 94, in
capped_gvs = [(tf.clip_by_value(grad, -1., 1.), var) for grad, var in gvs]
File "run.py", line 94, in
capped_gvs = [(tf.clip_by_value(grad, -1., 1.), var) for grad, var in gvs]
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\util\dispatch.py", line 180, in wrapper
return target(*args, **kwargs)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\ops\clip_ops.py", line 69, in clip_by_value
t = ops.convert_to_tensor(t, name="t")
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\ops.py", line 1039, in convert_to_tensor
return convert_to_tensor_v2(value, dtype, preferred_dtype, name)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\ops.py", line 1097, in convert_to_tensor_v2
as_ref=False)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\ops.py", line 1175, in internal_convert_to_tensor
ret = conversion_func(value, dtype=dtype, name=name, as_ref=as_ref)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\constant_op.py", line 304, in _constant_tensor_conversion_function
return constant(v, dtype=dtype, name=name)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\constant_op.py", line 245, in constant
allow_broadcast=True)
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\constant_op.py", line 283, in _constant_impl
allow_broadcast=allow_broadcast))
File "C:\Users\Alou\anaconda3\envs\research_gpu\lib\site-packages\tensorflow\python\framework\tensor_util.py", line 454, in make_tensor_proto
raise ValueError("None values not supported.")
ValueError: None values not supported.

Compatible with Tensorflow 2.0

Some of the functions, called in the tstDemo.py-file, e.g. reset_default_graph(), ConfigProto(), global_variables(), do not work for Tensorflow 2.0. A way around it is described here: https://www.tensorflow.org/api_docs/python/tf/compat/v1/reset_default_graph

Why do you calculate the gradient for the conjugate gradient mannually?

Hi, Hemant,

I'm diving into the implementation detail of conjugate gradient.
I was wondering why do you calculate the gradient for the conjugate gradient mannually instead of using auto diff by tensorflow ?

Thank you in advance!

parameter setting to reproduce the result in tstDemo.py

Hi, I have set the parameter as your parer write: layers N =5, iterations K =10 , others are the same in trn.py:
nLayers=5
epochs=100
batchSize=4
gradientMethod='MG'
K=10
sigma=0.01
both gradient method 'MG' and 'AG' have been tried, and variable reuse=tf.AUTO_REUSE also has been set. However, I can't reproduce the PSNR as high as 39dB in the provided model. My result is 33dB or less.

So, I want to know, how could I set the parameters in your code properly to reproduce your result?
Best !

ResourceExhaustedError: OOM when allocating tensor with shape[12,256,232]

Hi, I have several questions when training using trn.py
1.Out of memory
tensorflow.python.framework.errors_impl.ResourceExhaustedError: OOM when allocating tensor with shape[12,256,232]
My system is windows 10, 16 GB RAM, TIAN XP(12 GB). The model in your paper is relatively small and the dataset is also not big. I want to know if the dataset was readed into GPU once?

2.How to realize variable share in different iterations in your code:
`def createLayer(x, szW, trainning,lastLayer):

W=tf.get_variable('W',shape=szW,initializer=tf.contrib.layers.xavier_initializer())
x = tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME')
xbn=tf.layers.batch_normalization(x,training=trainning,fused=True,name='BN')

if not(lastLayer):
    return tf.nn.relu(xbn)
else:
    return xbn`

3.I want to use the two-step approach to initialize the network for training. The firstly trained model was different with the model with 10 iterations, using the following approach is not feasible
loadFirstTrainedChkPoint = tf.train.latest_checkpoint(FirstTrainedmodel) new_saver = tf.train.import_meta_graph(modelDir+'/model.meta') new_saver.restore(sess, loadFirstTrainedChkPoint)

Thanks!

Is the version of tensorflow above 1.7 or below?

I can run with tenflow-gpu==1.13.1, but the code can't end normally after the picture appears, that is, "Process finished with exit code 0" is not displayed, and using version 1.3 or 1.2 again will report an error

The knee data seems to be lost

Hi, Hemant,

The knee data in the link you provide in the Readme.md is invalid. Could you please check the data?

Thanks a lot !

The knee data you used for training

Hi!

  I saw you have added the demo of  knee data from the NYU, what a good job! 
  But there are many kinds of  knee data corresponding to different sequences in the link you provided, would you to give a introduction about which sequence you choose for training and what is the config for training. 
   By the way, the knee data of one sequence  contains about 20 persons and each person contains  20 or 30 image slices.  So, is there the necessary to use a per-trained K=1 model for initializing your network when reproducing the result provided in the  knee_demo_code.py?

Best!
TIeyuan

Question in Conjugate Gradient method

Hello Sir,

First of all, thank you for maintaining and sharing this repository. The idea of using the CG method for generating MR images is truly impressive 😃 .

My query was related to the Conjugate Gradient methods. In the file model.py, under the myCG function, at line 101,

alpha = rTr / tf.to_float(tf.reduce_sum(tf.conj(p)*Ap))

The tf.to_float() command discards the imaginary part. I checked the values over multiple steps and found that the imaginary values were very close to 0, so they could be discarded. It makes sense to discard the imaginary values when something like x*tf.conj(x) happens, as it would be = modulus(x)², which happens in line 105 when rtrNew is calculated. I'm just curious to know how did you reach the conclusion that these values would be small enough to be discarded 🟢 . Did you also check the values or was it based on some matrix properties of Ap? I read this topic in 'Numerical recipes for C' and it mentions that A should be symmetric and positive definite. But I don't know what properties Ap has which makes it behave the way it does in the code.

I would highly appreciate your feedback 😄