Final_Complex_rRAKI.py

import tensorflow as tf
tf.compat.v1.disable_eager_execution()
import scipy.io as sio
import numpy as np
import time
import os
import math
import matplotlib.pyplot as plt
os.environ["CUDA_VISIBLE_DEVICES"] = "0"

#KNEE FILES

# filename = 'file1000254'
# matfn = '/content/drive/My Drive/ICASSP_rRAKI/data/knee_multicoil_val/multicoil_val/' + filename + '.mat'
# name_weight = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Knee/' + filename + '_nocrop_weight.mat'
# recon_name = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Knee/' + filename + '_nocrop_recon.mat'
# loss_name = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Knee/' + filename + '_nocrop_loss.mat'
# tensorboard_dir = './graphs' + filename

#BRAIN FILES

filename = 'file_brain_AXFLAIR_206_6000241'
matfn = '/content/drive/My Drive/ICASSP_rRAKI/data/brain_multicoil_train/brain_multicoil_train/' + filename + '.mat'
name_weight = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Brain/' + filename + '_nocrop_weight.mat'
recon_name = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Brain/' + filename + '_nocrop_recon.mat'
loss_name = '/content/drive/My Drive/ICASSP_rRAKI/Complex_rRAKI_Results/Latest_pplanerelu/Brain/' + filename + '_nocrop_loss.mat'
tensorboard_dir = './graphs' + filename + '_tensorboard'

# def complex_to_channels(image):            
#     """Convert data from complex to channels."""
#     image_out = tf.stack([tf.math.real(image), tf.math.imag(image)], axis=-1)
#     shape_out = tf.concat(
#         [tf.shape(input=image)[:-1], [image.shape[-1] * 2]], axis=0)
#     image_out = tf.reshape(image_out, shape_out)
#     return image_out

def weight_variable(shape,vari_name):                   
    initial = tf.random.truncated_normal(shape, stddev=0.1,dtype=tf.float32)
    return tf.Variable(initial,name = vari_name)

def init_Newrelu_coefficients():
    coeff = tf.Variable(np.random.normal()) 
    coeff1 = tf.Variable(np.random.normal()) 
    coeff2 = tf.Variable(np.random.normal())
    return coeff, coeff1, coeff2
    
hyp_den = 3
        
def Newrelu(x):
    s = int(x.shape[-1]) 
    a, b, c = init_Newrelu_coefficients()
    prod = tf.add(tf.add(tf.multiply(x[:,:,:,:s//2], a), c), tf.multiply(x[:,:,:,s//2:], b))
    final_prod = tf.concat([prod, prod], axis=-1)
    return tf.where(tf.math.greater_equal(final_prod,0), x, (a+b+c)/hyp_den*x), a, b, c
    
def Newrelu1(x, a, b, c):
    s = int(x.shape[-1]) 
    prod = tf.add(tf.add(tf.multiply(x[:,:,:,:s//2], a), c), tf.multiply(x[:,:,:,s//2:], b))
    final_prod = tf.concat([prod, prod], axis=-1)
    return tf.where(tf.math.greater_equal(final_prod,0), x, (a+b+c)/hyp_den*x)

def cconv2d_dilate(x, real_W, imag_W, dilate_rate):
    shape_x = x.shape  #(1, 640, 25, 32)
    s = int(0 if shape_x[-1] is None else shape_x[-1])
    real_x = x[:,:,:,:s//2]
    imag_x = x[:,:,:,s//2:]

    real_real = tf.nn.convolution(input=real_x, filters=real_W, padding='VALID',dilations = [1,dilate_rate])
    real_imag = tf.nn.convolution(input=real_x, filters=imag_W, padding='VALID',dilations = [1,dilate_rate])
    imag_real = tf.nn.convolution(input=imag_x, filters=real_W, padding='VALID',dilations = [1,dilate_rate])
    imag_imag = tf.nn.convolution(input=imag_x, filters=imag_W, padding='VALID',dilations = [1,dilate_rate])   

    out_real = real_real-imag_imag
    out_imag = imag_real+real_imag

    # complex_output = tf.complex(out_real, out_imag)
    # channels_output = complex_to_channels(complex_output)
    
    channels_output = tf.concat([out_real, out_imag], axis=-1)

    return channels_output

def learning(accrate_input,sess,graph_dir,chnl):
    
    print("Final learning version ")
    
    input_ACS = tf.compat.v1.placeholder(tf.float32, [1, ACS_dim_X,ACS_dim_Y,ACS_dim_Z])                                   
    input_Target = tf.compat.v1.placeholder(tf.float32, [1, target_dim_X,target_dim_Y,target_dim_Z])           
    
    Input = tf.reshape(input_ACS, [1, ACS_dim_X, ACS_dim_Y, ACS_dim_Z])         

    [target_dim0,target_dim1,target_dim2,target_dim3] = np.shape(target)

    ker_conv_r = weight_variable([kernel_x_1, kernel_y_1, ACS_dim_Z//2, target_dim3//2],'G1')
    ker_conv_i = weight_variable([kernel_x_1, kernel_y_1, ACS_dim_Z//2, target_dim3//2],'G2')
    grp_conv = cconv2d_dilate(Input, ker_conv_r, ker_conv_i ,accrate_input)

    x_shift = np.int32(np.floor(kernel_last_x/2))
    
    [aa,bb,dim_yy,cc]=np.shape(grp_conv);

    grap_y_start = np.int32(  (np.ceil(kernel_y_2/2)-1) + (np.ceil(kernel_last_y/2)-1)) * accrate_input
    grap_y_end = np.int32(dim_yy) - np.int32(( (np.floor(kernel_y_2/2)) + np.floor(kernel_last_y/2))) * accrate_input - 1
    grap_y_start = np.int32(grap_y_start);
    grap_y_end = np.int32(grap_y_end+1);

    grapRes =  grp_conv[:,x_shift:x_shift+target_dim_X,grap_y_start:grap_y_end,:];
    
    #conv1 layer 
    W_conv1_r = weight_variable([kernel_x_1, kernel_y_1, ACS_dim_Z//2, layer1_channels//2],'W11')
    W_conv1_i = weight_variable([kernel_x_1, kernel_y_1, ACS_dim_Z//2, layer1_channels//2],'W12') 
    out1 = Newrelu(cconv2d_dilate(Input, W_conv1_r, W_conv1_i, accrate_input))
    h_conv1 = out1[0]
    coeff_a1 = out1[1]
    coeff_b1 = out1[2]
    coeff_c1 = out1[3]

    # conv2 layer
    W_conv2_r = weight_variable([kernel_x_2, kernel_y_2, layer1_channels//2, layer2_channels//2],'W21')
    W_conv2_i = weight_variable([kernel_x_2, kernel_y_2, layer1_channels//2, layer2_channels//2],'W22')
    out2 = Newrelu(cconv2d_dilate(h_conv1, W_conv2_r, W_conv2_i, accrate_input))
    h_conv2 = out2[0]
    coeff_a2 = out2[1]
    coeff_b2 = out2[2]
    coeff_c2 = out2[3]

    ## conv3 layer
    W_conv3_r = weight_variable([kernel_last_x, kernel_last_y, layer2_channels//2, target_dim3//2],'W31')
    W_conv3_i = weight_variable([kernel_last_x, kernel_last_y, layer2_channels//2, target_dim3//2],'W32')
    h_conv3 = cconv2d_dilate(h_conv2, W_conv3_r, W_conv3_i, accrate_input)

    error_norm = 1*tf.norm(tensor=input_Target - grapRes - h_conv3) + 1*tf.norm(tensor=input_Target - grapRes) 
    
    print("MOMENTUM LEARNING")
    train_step = tf.compat.v1.train.MomentumOptimizer(learning_rate=0.001,momentum=0.9,use_nesterov=False).minimize(error_norm)
    
    if int((tf.__version__).split('.')[1]) < 12 and int((tf.__version__).split('.')[0]) < 1:
        init = tf.compat.v1.initialize_all_variables()
    else:
        init = tf.compat.v1.global_variables_initializer()
    sess.run(init)
    
    print("INIT act coefficients ", sess.run(coeff_a1),sess.run(coeff_b1),sess.run(coeff_c1), sess.run(coeff_a2),sess.run(coeff_b2),sess.run(coeff_c2))

    # writer = tf.compat.v1.summary.FileWriter(graph_dir)
    writer = tf.summary.create_file_writer(graph_dir)
    
    with writer.as_default():
        for i in range(epos):
            
            sess.run(train_step, feed_dict={input_ACS: ACS, input_Target: target})
            # summary = tf.compat.v1.Summary(value=[tf.compat.v1.Summary.Value(value=sess.run(error_norm,feed_dict={input_ACS: ACS, input_Target: target}))])
            # tf. summary.scalar("error_norm", sess.run(error_norm,feed_dict={input_ACS: ACS, input_Target: target}), step=i)
            # writer.add_summary(summary, i)
            if i%10 == 0:
                index = int(i//10)
                error_arr[chnl, index] = sess.run(error_norm,feed_dict={input_ACS: ACS, input_Target: target})
            if i % 100 == 0:                                                                         
                error_now=sess.run(error_norm,feed_dict={input_ACS: ACS, input_Target: target})     
                print('The',i,'th iteration gives an error',error_now)  
            
    writer.flush() # make sure everything is written to disk
    # writer.close()
            
    error = sess.run(error_norm,feed_dict={input_ACS: ACS, input_Target: target})
    print("LEARNT act coefficients ", sess.run(coeff_a1),sess.run(coeff_b1),sess.run(coeff_c1), sess.run(coeff_a2),sess.run(coeff_b2),sess.run(coeff_c2))
    return [sess.run(ker_conv_r),sess.run(ker_conv_i),sess.run(W_conv1_r),sess.run(W_conv1_i),sess.run(W_conv2_r),sess.run(W_conv2_i),sess.run(W_conv3_r),sess.run(W_conv3_i),sess.run(coeff_a1),sess.run(coeff_b1),sess.run(coeff_c1),sess.run(coeff_a2),sess.run(coeff_b2),sess.run(coeff_c2),error]   
  
                                                            
def cnn_3layer(input_kspace,gkerr,gkeri,w1r,w1i,w2r,w2i,w3r,w3i,acc_rate,cc1,cc2,sess):                 
    grap = cconv2d_dilate(input_kspace, gkerr, gkeri ,acc_rate)
    x_shift = np.int32(np.floor(kernel_last_x/2))
    
    [aa,dim_x,dim_yy,dd] = np.shape(grap);

    grap_y_start = np.int32((np.ceil(kernel_y_2/2)-1) + (np.ceil(kernel_last_y/2)-1)) * acc_rate
    grap_y_end = np.int32(dim_yy) - np.int32(( (np.floor(kernel_y_2/2)) + np.floor(kernel_last_y/2))) * acc_rate - 1

    grap_y_start = np.int32(grap_y_start);
    grap_y_end = np.int32(grap_y_end+1);

    effectiveGrappa =  grap[:, x_shift:dim_x-x_shift, grap_y_start:grap_y_end, :];
    
    cc1 = np.squeeze(cc1)
    cc2 = np.squeeze(cc2)
    
    print("Coefficients to layer1 in reconstruction ", cc1[0], cc1[1], cc1[2])
    h_conv1 = Newrelu1(cconv2d_dilate(input_kspace, w1r, w1i, acc_rate), cc1[0], cc1[1], cc1[2]) 
    print("Coefficients to layer2 in reconstruction ", cc2[0], cc2[1], cc2[2])
    h_conv2 = Newrelu1(cconv2d_dilate(h_conv1, w2r, w2i, acc_rate), cc2[0], cc2[1], cc2[2])
    h_conv3 = cconv2d_dilate(h_conv2, w3r, w3i, acc_rate) 
    
    print('grap res shape = ',np.shape(grap))
    print('effective grap shape = ',np.shape(effectiveGrappa))
    print('h_conv shape = ',np.shape(h_conv3))
    return sess.run(effectiveGrappa+h_conv3), sess.run(effectiveGrappa),sess.run(h_conv3), sess.run(grap) 


#Initialization


no_ACS_flag = 0
rate = 5
ACSrange = 24
phaseshiftflag = 0;
kspace = sio.loadmat(matfn)
kspace = kspace['kspace'] 
mask = np.zeros_like(kspace)
[row,col,coil] = mask.shape
mask[:,::rate,:] = 1
midpoi = col//2
mask[:,midpoi-ACSrange//2+1:midpoi+ACSrange//2,:]=1
kspace = kspace*mask
normalize = 1/np.max(abs(kspace[:]))
kspace = np.multiply(kspace,normalize)    

[m1,n1,no_ch] = np.shape(kspace) 
no_inds = 1

kspace_all = kspace;
kx = np.transpose(np.int32([(range(1,m1+1))]))                          
ky = np.int32([(range(1,n1+1))])

if phaseshiftflag ==1:
    phase_shifts = np.dot(np.exp(-1j * 2 * 3.1415926535 / m1 * (m1/2-1) * kx ),np.exp(-1j * 2 * 3.14159265358979 / n1 * (n1/2-1) * ky ))
    for channel in range(0,no_ch):
        kspace_all[:,:,channel] = np.multiply(kspace_all[:,:,channel],phase_shifts)

kspace = np.copy(kspace_all)
mask = np.squeeze(np.sum(np.sum(np.abs(kspace),axis=0),axis=1))>0
picks = np.where(mask == 1);                                  

kspace_NEVER_TOUCH = np.copy(kspace_all)

mask = np.squeeze(np.sum(np.sum(np.abs(kspace),axis=0),axis=1))>0;  
picks = np.where(mask == 1);                                  
d_picks = np.diff(picks,1)  
indic = np.where(d_picks == 1);

mask_x = np.squeeze(np.sum(np.sum(np.abs(kspace),axis=2),axis=1))>0;
picks_x = np.where(mask_x == 1);
x_start = picks_x[0][0]
x_end = picks_x[0][-1]

print('ACS signal found in the input data')
indic = indic[1][:]
center_start = picks[0][indic[0]];
center_end = picks[0][indic[-1]+1];

print('START AND END OF ACS ',center_start,center_end)

ACS = kspace[x_start:x_end+1,center_start:center_end+1,:]
[ACS_dim_X, ACS_dim_Y, ACS_dim_Z] = np.shape(ACS)
ACS_re = np.zeros([ACS_dim_X,ACS_dim_Y,ACS_dim_Z*2])
ACS_re[:,:,0:no_ch] = np.real(ACS)
ACS_re[:,:,no_ch:no_ch*2] = np.imag(ACS)


acc_rate = d_picks[0][0]
no_channels = ACS_dim_Z*2


kernel_x_1 = 5
kernel_y_1 = 2

kernel_x_2 = 1
kernel_y_2 = 1

kernel_last_x = 3
kernel_last_y = 2

layer1_channels = 32
layer2_channels = 8

print('layer1_channels, layer2_channels', layer1_channels, layer2_channels)


gker_all_r = np.zeros([kernel_x_1, kernel_y_1, no_channels//2, acc_rate - 1, no_channels//2],dtype=np.float32)
w1_all_r = np.zeros([kernel_x_1, kernel_y_1, no_channels//2, layer1_channels//2, no_channels//2],dtype=np.float32)
w2_all_r = np.zeros([kernel_x_2, kernel_y_2, layer1_channels//2,layer2_channels//2,no_channels//2],dtype=np.float32)
w3_all_r = np.zeros([kernel_last_x, kernel_last_y, layer2_channels//2,acc_rate - 1, no_channels//2],dtype=np.float32)   
coeffs1 = np.zeros([no_channels//2, 3], dtype=np.float32)

gker_all_i = np.zeros([kernel_x_1, kernel_y_1, no_channels//2, acc_rate - 1, no_channels//2],dtype=np.float32)
w1_all_i = np.zeros([kernel_x_1, kernel_y_1, no_channels//2, layer1_channels//2, no_channels//2],dtype=np.float32)
w2_all_i = np.zeros([kernel_x_2, kernel_y_2, layer1_channels//2,layer2_channels//2,no_channels//2],dtype=np.float32)
w3_all_i = np.zeros([kernel_last_x, kernel_last_y, layer2_channels//2, acc_rate - 1, no_channels//2],dtype=np.float32)
coeffs2 = np.zeros([no_channels//2, 3], dtype=np.float32)


target_x_start = np.int32(np.ceil(kernel_x_1/2) + np.floor(kernel_x_2/2) + np.floor(kernel_last_x/2) -1);
target_x_end = np.int32(ACS_dim_X - target_x_start -1)

#Learning

time_ALL_start = time.time()

[ACS_dim_X, ACS_dim_Y, ACS_dim_Z] = np.shape(ACS_re)
ACS = np.reshape(ACS_re, [1,ACS_dim_X, ACS_dim_Y, ACS_dim_Z]) # Batch, X, Y, Z
ACS = np.float32(ACS)  # here we use ACS instead of ACS_re for convenience

target_y_start = np.int32((np.ceil(kernel_y_1/2)-1) + (np.ceil(kernel_y_2/2)-1) + (np.ceil(kernel_last_y/2)-1)) * acc_rate;     
target_y_end = ACS_dim_Y  - np.int32((np.floor(kernel_y_1/2) + np.floor(kernel_y_2/2) + np.floor(kernel_last_y/2))) * acc_rate -1;

target_dim_X = target_x_end - target_x_start + 1
target_dim_Y = target_y_end - target_y_start + 1
target_dim_Z = (acc_rate - 1)*2 # *2 for complex valued 

print('go!')
time_Learn_start = time.time() # set timer

errorSum = 0;
config = tf.compat.v1.ConfigProto()
config.gpu_options.per_process_gpu_memory_fraction = 1/3; # avoid fully allocating.

epos = 2000
error_arr = np.zeros([ACS_dim_Z//2,epos//10])

for ind_c in range(ACS_dim_Z//2):
    sess = tf.compat.v1.Session(config=config)  
    target = np.zeros([1,target_dim_X,target_dim_Y,target_dim_Z])
    print('learning channel #',ind_c+1)
    time_channel_start = time.time()
    lim = acc_rate-1 
    for ind_acc in range(lim):
        target_y_start = np.int32((np.ceil(kernel_y_1/2)-1) + (np.ceil(kernel_y_2/2)-1) + (np.ceil(kernel_last_y/2)-1)) * acc_rate + ind_acc + 1
        target_y_end = ACS_dim_Y  - np.int32((np.floor(kernel_y_1/2) + (np.floor(kernel_y_2/2)) + np.floor(kernel_last_y/2))) * acc_rate + ind_acc
        target[0,:,:,ind_acc] = ACS[0,target_x_start:target_x_end + 1, target_y_start:target_y_end +1,ind_c];
        target[0,:,:,ind_acc+lim] = ACS[0,target_x_start:target_x_end + 1, target_y_start:target_y_end +1,ind_c+no_channels//2];

    graph_dir = tensorboard_dir + str(ind_c)
    [gkerr,gkeri,w1r,w1i,w2r,w2i,w3r,w3i,error,coeff_a1l,coeff_b1l,coeff_c1l,coeff_a2l,coeff_b2l,coeff_c2l]=learning(acc_rate,sess,graph_dir,ind_c) 
    gker_all_r[:,:,:,:,ind_c]=gkerr
    w1_all_r[:,:,:,:,ind_c] = w1r
    w2_all_r[:,:,:,:,ind_c] = w2r
    w3_all_r[:,:,:,:,ind_c] = w3r
    coeffs1[ind_c, 0] = coeff_a1l
    coeffs1[ind_c, 1] = coeff_b1l
    coeffs1[ind_c, 2] = coeff_c1l

    gker_all_i[:,:,:,:,ind_c]=gkeri
    w1_all_i[:,:,:,:,ind_c] = w1i
    w2_all_i[:,:,:,:,ind_c] = w2i
    w3_all_i[:,:,:,:,ind_c] = w3i 
    coeffs2[ind_c, 0] = coeff_a2l
    coeffs2[ind_c, 1] = coeff_b2l
    coeffs2[ind_c, 2] = coeff_c2l
    
    time_channel_end = time.time()
    print('Time Cost:',time_channel_end-time_channel_start,'s')
    # print('Norm of Error = ',error)
    errorSum = errorSum + error

    sess.close()
    tf.compat.v1.reset_default_graph()

time_Learn_end = time.time();
print('lerning step costs:',(time_Learn_end - time_Learn_start)/60,'min') # get time

sio.savemat(name_weight, {'gkerr':gker_all_r,'w1r': w1_all_r,'w2r': w2_all_r,'w3r': w3_all_r, 'gkeri':gker_all_i,'w1i': w1_all_i,'w2i': w2_all_i,'w3i': w3_all_i, 'coeffs1' : coeffs1, 'coeffs2' : coeffs2})  
np.save(loss_name, error_arr)


weightfile = sio.loadmat(name_weight)
gker_allr = weightfile['gkerr'] # get kspace
gker_alli = weightfile['gkeri'] # get kspace
w1_allr = weightfile['w1r'] # get kspace
w1_alli = weightfile['w1i'] # get kspace
w2_allr = weightfile['w2r'] # get kspace
w2_alli = weightfile['w2i'] # get kspace
w3_allr = weightfile['w3r'] # get kspace
w3_alli = weightfile['w3i'] # get kspace
c1 = weightfile['coeffs1']
c2 = weightfile['coeffs2']
print("Shapes of C1, C2 ", c1.shape, c2.shape)

#Reconstruction

kspace_recon_all = np.copy(kspace_all)
kspace_recon_all_nocenter = np.copy(kspace_all)

kspace = np.copy(kspace_all)

over_samp = np.setdiff1d(picks,np.int32([range(0, n1,acc_rate)]))
kspace_und = kspace
kspace_und[:,over_samp,:] = 0;
[dim_kspaceUnd_X,dim_kspaceUnd_Y,dim_kspaceUnd_Z] = np.shape(kspace_und)

kspace_und_re = np.zeros([dim_kspaceUnd_X,dim_kspaceUnd_Y,dim_kspaceUnd_Z*2])
kspace_und_re[:,:,0:dim_kspaceUnd_Z] = np.real(kspace_und)
kspace_und_re[:,:,dim_kspaceUnd_Z:dim_kspaceUnd_Z*2] = np.imag(kspace_und)
kspace_und_re = np.float32(kspace_und_re)
kspace_und_re = np.reshape(kspace_und_re,[1,dim_kspaceUnd_X,dim_kspaceUnd_Y,dim_kspaceUnd_Z*2])
kspace_recon = kspace_und_re

raki_recon = np.zeros_like(kspace_recon)
grap_recon = np.copy(kspace_recon)

config = tf.compat.v1.ConfigProto()
config.gpu_options.per_process_gpu_memory_fraction = 1/3 ; 

for ind_c in range(0,no_channels//2):
    print('Reconstructing Channel #',ind_c+1)


    sess = tf.compat.v1.Session(config=config)  # tensorflow initialize
    if int((tf.__version__).split('.')[1]) < 12 and int((tf.__version__).split('.')[0]) < 1:
        init = tf.compat.v1.initialize_all_variables()
    else:
        init = tf.compat.v1.global_variables_initializer()
    sess.run(init)

    gkerr = np.float32(gker_allr[:,:,:,:,ind_c])
    gkeri = np.float32(gker_alli[:,:,:,:,ind_c])
    w1r = np.float32(w1_allr[:,:,:,:,ind_c])
    w1i = np.float32(w1_alli[:,:,:,:,ind_c])
    w2r = np.float32(w2_allr[:,:,:,:,ind_c])
    w2i = np.float32(w2_alli[:,:,:,:,ind_c])
    w3r = np.float32(w3_allr[:,:,:,:,ind_c])
    w3i = np.float32(w3_alli[:,:,:,:,ind_c])
    curr_c1 = np.float32(c1[ind_c,:])
    curr_c2 = np.float32(c2[ind_c,:])               
        
    [res,grap,raki,rawgrap] = cnn_3layer(kspace_und_re,gkerr,gkeri,w1r,w1i,w2r,w2i,w3r,w3i,acc_rate,curr_c1,curr_c2,sess)  
    target_x_end_kspace = dim_kspaceUnd_X - target_x_start;
    
    for ind_acc in range(0,acc_rate-1):

        target_y_start = np.int32((np.ceil(kernel_y_1/2)-1) + np.int32((np.ceil(kernel_y_2/2)-1)) + np.int32(np.ceil(kernel_last_y/2)-1)) * acc_rate + ind_acc + 1
        target_y_end_kspace = dim_kspaceUnd_Y - np.int32((np.floor(kernel_y_1/2)) + (np.floor(kernel_y_2/2)) + np.floor(kernel_last_y/2)) * acc_rate + ind_acc;
        kspace_recon[0,target_x_start:target_x_end_kspace,target_y_start:target_y_end_kspace+1:acc_rate,ind_c] = res[0,:,::acc_rate,ind_acc]

        grap_recon[0,target_x_start:target_x_end_kspace,target_y_start:target_y_end_kspace+1:acc_rate,ind_c] = grap[0,:,::acc_rate,ind_acc]

        raki_recon[0,target_x_start:target_x_end_kspace,target_y_start:target_y_end_kspace+1:acc_rate,ind_c] = raki[0,:,::acc_rate,ind_acc]


    print('total parameters ', np.sum([np.prod(v.get_shape().as_list()) for v in tf.compat.v1.trainable_variables()]))
    sess.close()
    tf.compat.v1.reset_default_graph()

kspace_recon = np.squeeze(kspace_recon)

kspace_recon_complex = (kspace_recon[:,:,0:np.int32(no_channels/2)] + np.multiply(kspace_recon[:,:,np.int32(no_channels/2):no_channels],1j))
kspace_recon_all_nocenter[:,:,:] = np.copy(kspace_recon_complex); # im_ind = 1, skip one dim

grap_recon = np.squeeze(grap_recon)

grap_recon_complex = (grap_recon[:,:,0:np.int32(no_channels/2)] + np.multiply(grap_recon[:,:,np.int32(no_channels/2):no_channels],1j))

raki_recon = np.squeeze(raki_recon)

raki_recon_complex = (raki_recon[:,:,0:np.int32(no_channels/2)] + np.multiply(raki_recon[:,:,np.int32(no_channels/2):no_channels],1j)) 

if no_ACS_flag == 0:  # if we have ACS region in kspace, put them back
    kspace_recon_complex[:,center_start:center_end,:] = kspace_NEVER_TOUCH[:,center_start:center_end,:]
    print('ACS signal has been put back')
else:
    print('No ACS signal is put into k-space')

kspace_recon_all[:,:,:] = kspace_recon_complex; # im_ind = 1, skip one dim

for sli in range(0,no_ch):
    kspace_recon_all[:,:,sli] = np.fft.ifft2(kspace_recon_all[:,:,sli])

rssq = (np.sum(np.abs(kspace_recon_all)**2,2)**(0.5))
sio.savemat(recon_name,{'kspace_all':kspace_recon_complex,'kspace_all_noACS':kspace_recon_all_nocenter,'grap_all':grap_recon_complex,'raki_all':raki_recon_complex,'rawgrap':rawgrap,'effectiveGrap':grap})  # save the results

time_ALL_end = time.time()
print('All process costs ',(time_ALL_end-time_ALL_start)/60,'mins')
print('Error Average in Training is ',errorSum/no_channels)