<jupyter_start><jupyter_text>This is an MRI based reconstruction demo, for 2D MRI data. The network is relatively similar to the AutoMap technique (https://arxiv.org/abs/1704.08841). This is a relatively 'brute force' aproach to image reconstruction in which the transoform is given no direct knowledge of the physics (although the network architecture is a bit tuned to the problem). In this work, we are assuming one direction is fully sampled (i.e. frequency encoded). MRI SamplingIn MRI the data is often discretely Fourier transoformed in one direction leading to the discretized signal model:$s(k)=\sum_{j=1}^{N}\rho (x_j)e^{i2\pi kx}$The expected reconstruction for fully sampled data is an inverse discrete Fourier transform: $s(x)=\sum_{j=1}^{N}s(k_j)e^{i2\pi k_j x}$ Questions to think about:1) What is the minimal network architecture to compute a DFT? It's a square matrix multiply.2) What is the apropriate loss function if we wish to train an image reconstruction? 3) What is the role of the convolutional layers? When are they needed?4) What is the network learning if you train on simulated images? PyTorchThis exercise is written in PyTorch (https://pytorch.org/docs/stable/index.html). There is also a Keras version available online. PyTorch has some different coding constructs but there are similar ideas to Keras.   Setup<jupyter_code>''' 
In python you need to import libraries in order to use them. 
'''

# Using PyTorch for this code
import torch
  
# Utilities
import numpy as np
import math 

# Plotting
import matplotlib.pyplot as plt
from IPython.display import clear_output

''' 
We will define some functions to use later. Normally you might put these in another file and import them
'''

# Some support functions
def montage( img_in, size=(3,5) ):
    for j in range(size[0]):
        plot_image = img_in[0+size[1]*j,:,:]
        for i in range(size[1]-1):
            plot_image = np.concatenate((plot_image, img_in[1+i+size[1]*j,:,:]), axis=1)
        
        if j == 0:
            img = plot_image
        else:
            img = np.concatenate((img,plot_image),axis=0)
    return img
  
def complex_to_channels( img_in):
  return(np.stack(img_in.real,img_in.imag))
  
def channels_to_complex( img_in):
  return(img_in[...,0]+1j*img_in[...,1])


'''
Mount your google drive, we'll grab data from the shared folder
'''
from google.colab import drive
drive.mount('/content/drive')<jupyter_output><empty_output><jupyter_text>Training ImagesWe are training this network using a simulation enviroment. Images are grabbed from a set of MRI brain images. This example is using the output dicoms which have been preprocessed. We then simulate the image to MRI raw data conversion.<jupyter_code># load training, validation, and testing data
import h5py
with h5py.File('/content/drive/MyDrive/ML4MI_BOOTCAMP_DATA/ImageReconstructionData.h5','r') as hf:
  x_train = np.array(hf['x_train'])
  x_val = np.array(hf['x_val'])

print(f'Validate Dataset Size {x_val.shape}')
print(f'Train Dataset Size {x_train.shape}')
N = x_train.shape[-2]<jupyter_output><empty_output><jupyter_text>Simulate SamplingMRI data generation is aproximately dsicrete sampling of a continous Fourier transform the the data. In this example, we are using a Discrete Fourier transform to aproximate this. We also consider the case when we randomly remove data points. This would allow us to go faster and is used in compressed sensing application ( e.g. https://onlinelibrary.wiley.com/doi/pdf/10.1002/mrm.21391 ). Noise is added a complex, white, gaussian noise (MRI noise is so called Johnson/Nyquist noise). Things to try:1) Add higher levels of noise. What happens to the training rate and output images? 2) Increase the undersampling rate. How does the neural network compare to traditional aproaches? 3) Comment the FFT shift, does the network still learn the transform?<jupyter_code>'''
The creates a sampling mask which can be used to subsample the data.
'''

# Get the number of phase encodes
undersample_factor = 1.5 
noise_level = 0.001

number_phase_encodes = int(N/undersample_factor)
print('Using ' + str(number_phase_encodes) + ' phase encode')

# Create a random mask to resample the data
idx = np.full(N, False)
idx[:number_phase_encodes] = True
np.random.seed(1) # Keep this one so code is reproducible
np.random.shuffle(idx)
sampling_mask = idx


'''
Fourier transform, subsample, and add noise (Train Data)
'''
Nexamples = x_train.shape[0]
kspace_train = np.zeros((Nexamples,N,number_phase_encodes,2),x_train.dtype)
for example in range(x_train.shape[0]):
  
  if example % 1000 == 0:
    print(f'Working on example {example} of {Nexamples}')
      
  # Grab one image
  temp = x_train[example,:,:,0] + 1j*x_train[example,:,:,1]
  
  # Fourier Transform
  kspace_temp = np.fft.fftn(temp,axes=(1,))/N
  kspace_temp = np.fft.fftshift(kspace_temp,axes=(1,))
  kspace_temp =np.stack( (kspace_temp.real, kspace_temp.imag), axis=-1)
  
  # Subsample
  kspace_temp = kspace_temp[:,sampling_mask,:]
  
  # Add noise
  kspace_temp += noise_level*np.random.randn(*kspace_temp.shape)

  # Put back
  kspace_train[example,:,:,:] = kspace_temp
  
print('Dimensions of training data are ' + str(kspace_train.shape) + '[ Examples x Nx x Ny x Channels]')


'''
Fourier transform, subsample, and add noise (Validate Data) 
'''
Nexamples = x_val.shape[0]
kspace_val = np.zeros((Nexamples,N,number_phase_encodes,2),x_train.dtype)
for example in range(x_val.shape[0]):
  
  if example % 1000 == 0:
    print(f'Working on example {example} of {Nexamples}')
      
  # Grab one image
  temp = x_val[example,:,:,0] + 1j*x_val[example,:,:,1]
  
  # Fourier Transform
  kspace_temp = np.fft.fftn(temp,axes=(1,))/N
  kspace_temp = np.fft.fftshift(kspace_temp,axes=(1,))
  kspace_temp =np.stack( (kspace_temp.real, kspace_temp.imag), axis=-1)
  
  # Subsample
  kspace_temp = kspace_temp[:,sampling_mask,:]
  
  # Add noise
  kspace_temp += noise_level*np.random.randn(*kspace_temp.shape)

  # Put back
  kspace_val[example,:,:,:] = kspace_temp
  
print('Dimensions of validation data are ' + str(kspace_val.shape) + '[ Examples x Nx x Ny x Channels]')
  
example = 400

# Show one image and k-space pair
img = x_train[example,:,:,0] + 1j*x_train[example,:,:,1]

plt.figure()
plt.subplot(121)
plt.imshow(np.abs(img),cmap='gray')
plt.title('Grayscale')

img = kspace_train[example,:,:,0] + 1j*kspace_train[example,:,:,1]

plt.subplot(122)
plt.imshow(np.abs(img),cmap='gray')
plt.title('K-Space (1D FFT)')
plt.show()<jupyter_output><empty_output><jupyter_text>Build the network architectureIn PyTorch, we need to create a model just as we do in Keras. However, models are usually defined as a class. A class is a structure that holds data and functions. We base our classes on torch.nn.Module and will only modify the __init__ and forward functions. The __init__ module will define the modules (convultions, activations) that we will use in the forward pass. This method is very flexible and PyTorch's autograd functionally will generally handle all the backpropagation. **FYI** PyTorch defaults to channels first. So the images will be [batch index, channels, x, y]We will define the following functions.*   FullyConnected a modules to perform the fully connected layers*   Denoiser a module to perform denoising using convolutions*   AutoMap a simple module which sequentially applies FullyConnected and Denoiser<jupyter_code>import torch
import torch.nn as nn
import torch.nn.functional as F

class FullyConnected(torch.nn.Module):
  """A fully connected network to transform from k-space to images.
  """

  def __init__(self, input_phase_encodes, image_size):
    super(FullyConnected, self).__init__()
    self.input_phase_encodes = input_phase_encodes
    self.image_size = image_size

    # These two convolution will be fully connected along one dimension. We are not
    # using a fully connected network in 2D due to memory contraints and the fact 
    # that the data is undersampling in one dimension. 
    self.conv1 = nn.Conv2d(in_channels=2, out_channels=2*image_size, kernel_size=(1, input_phase_encodes), padding='valid')
    self.conv2 = nn.Conv2d(in_channels=2, out_channels=2*image_size, kernel_size=(1, image_size), padding='valid')

  def forward(self, data):

    # This is a fully connected network followed by reshaping into an image
    layer = self.conv1(data)
    layer = layer.view((-1,2,self.image_size,self.image_size))
    layer = torch.moveaxis(layer,-2,-1)

    # This is a fully connected network followed by reshaping into an image
    layer = self.conv2(layer)
    layer = layer.view((-1,2,self.image_size,self.image_size))
    layer = torch.moveaxis(layer,-2,-1)

    return layer

class Denoiser(torch.nn.Module):
  """An image to image denoiser
  """

  def __init__(self, depth=3, initial_features=8):
    super(Denoiser, self).__init__()

    # The channels start as complex (2 channel) and we can increase them during the denoiser
    in_channels = 2
    out_channels = initial_features

    # This will build the list of operators, similar to keras sequential which each being a convolution followed by activation
    layers = []
    for l in range(depth):
      layers.append( nn.Conv2d(in_channels=in_channels, out_channels=out_channels, kernel_size=(3,3), padding='same'))
      in_channels = out_channels
      layers.append(torch.nn.ReLU())

    # Last layer is a convolution without activiation
    layers.append( nn.Conv2d(in_channels=in_channels, out_channels=2, kernel_size=(3,3), padding='same'))
    
    # This stores the layers in a list that will be tracked for backpropogation
    self.convolutions = nn.ModuleList(layers)
    
  def forward(self, image):
    for l in self.convolutions:
      image = l(image)

    return image

class AutoMap(torch.nn.Module):
  """A network combining the fully connected network with the image denoiser.
  """

  def __init__(self, input_phase_encodes, image_size, depth=3):
    super(AutoMap, self).__init__()

    self.fc_net = FullyConnected( input_phase_encodes, image_size)
    self.denoiser = Denoiser(depth)
    
  def forward(self, data):
    image = self.fc_net(data)
    image = self.denoiser(image)

    return image
   

# To run the model we need to create an object based on those class definitions above
recon_model = AutoMap(input_phase_encodes=number_phase_encodes, image_size=N)

# Models have a device they run on. We need to put the model on the gpu
recon_model = recon_model.cuda()

# Torch summary enables similar formatting to Keras
from torchsummary import summary
summary(recon_model, (2,120,120))<jupyter_output><empty_output><jupyter_text>Data LoaderIn PyTorch, it is very useful to have data as a torch Dataset. This could select all the images in a folder, generate data on the fly, etc. In this case, we will write a loader which just grabs data. Following this, we define a data loader which handles batching, shuffling, etc.<jupyter_code>class Dataset(torch.utils.data.Dataset):
  
  def __init__(self, kspace_data, image_data):
        self.kspace_data = np.moveaxis( kspace_data, -1, 1)
        self.image_data = np.moveaxis( image_data, -1, 1)

  def __len__(self):
        return self.kspace_data.shape[0]

  def __getitem__(self, idx):
        return self.kspace_data[idx], self.image_data[idx]

# Create datasets
dataset_train = Dataset( kspace_train, x_train)
dataset_val = Dataset( kspace_val, x_val)

# Create data loader to handle shuffling, batching, etc
train_generator = torch.utils.data.DataLoader(dataset_train, batch_size=32, shuffle=True)
val_generator = torch.utils.data.DataLoader(dataset_val, batch_size=32, shuffle=False)<jupyter_output><empty_output><jupyter_text>TrainingTraining in PyTorch is a bit more transparent than in Keras. The below will train the model with the key steps being:*   Define an optimizer which will optimize the paramaters we pass to it*   Define a loss function. This could be anything written in PyTorch but we will use a predefined function for mean squared error*   Loop over the data with a for loop.For each batch, we will:* Zero the gradients (gradients accumulate over multiple passes)* Perform a forward pass* Calculate the loss* Perform a backwards pass* Increment the weights using optimizer.step()<jupyter_code># Define an optimizer 
optimizer = torch.optim.Adam( recon_model.parameters(), lr=1e-3)

# Define a loss function
loss_fcn = nn.MSELoss()

# Get a device
use_cuda = torch.cuda.is_available()
device = torch.device("cuda:0" if use_cuda else "cpu")

# Ensure the model weights are on the desired device
recon_model.to(device)

# Define number of epochs
n_epochs = 20

# Empty list to store losses over epochs
train_losses = []
val_losses = []

# Epoch loop
for epoch in range(n_epochs):

  # Put the model in train mode
  recon_model.train()

  # Loop over training batches
  train_loss_avg = 0.0
  for kspace, image in train_generator:
    
    # Move data to the GPU
    image = image.to(device)
    kspace = kspace.to(device)

    # Zero out the gradients
    optimizer.zero_grad()

    # Forward pass
    image_guess = recon_model(kspace)

    # Loss
    loss = loss_fcn( image_guess, image)

    # Store loss 
    train_loss_avg += loss.detach()

    # Backwards with perform back propogation to get weights
    loss.backward()
    
    # Take a gradient descent step
    optimizer.step()

  # Store the average loss
  train_loss_avg /= len(train_generator)
  train_losses.append(train_loss_avg)

  # Switch to eval mode (weights fixed)
  recon_model.eval()

  # Loop eval batches
  val_loss_avg = 0.0
  for count, (kspace, image) in enumerate(val_generator):

    # Move data to the GPU
    image = image.to(device)
    kspace = kspace.to(device)
   
    # Forward pass
    image_guess = recon_model(kspace)

    # Loss
    loss = loss_fcn( image_guess, image)
    val_loss_avg += loss.detach()

    # Plotting images
    if count == 0:
      image_mag = torch.abs(image[0,0,...] + 1j*image[0,1,...])
      image_guess_mag = torch.abs(image_guess[0,0,...] + 1j*image_guess[0,1,...])

      clear_output(wait=True)
      plt.figure(figsize=(10,3))
      plt.subplot(131)
      plt.imshow(image_mag.detach().cpu().numpy(),cmap='gray')
      plt.title('Truth')
      plt.subplot(132)
      plt.imshow(image_guess_mag.detach().cpu().numpy(),cmap='gray')
      plt.title('NN Guess')

  # Store the average loss
  val_loss_avg /= len(val_generator)
  val_losses.append(val_loss_avg)

  # Plotting loss curve
  plt.subplot(133)
  plt.semilogy(train_losses, label="Loss")
  plt.semilogy(val_losses, label="Loss (validation)")
  plt.legend()
  plt.show()

  print(f'Epoch = {epoch} Loss = {train_loss_avg}, Val Loss = {val_loss_avg}')<jupyter_output><empty_output><jupyter_text>Run the model fit<jupyter_code>example = 400

# Test with synthetic data
kspace = kspace_val[example,...]
kspace = np.expand_dims(kspace,0)

image = x_val[example,...]
image = np.expand_dims(image,0)

# Get the prediction
kspace_torch = torch.tensor(np.moveaxis( kspace, -1, 1))
predicted_image = recon_model(kspace_torch.to(device))
predicted_image = predicted_image.detach().cpu().numpy()
predicted_image = np.moveaxis(predicted_image,1,-1)

# Convert to complex
predicted_image = np.squeeze(channels_to_complex(predicted_image))
act_image =  np.squeeze(channels_to_complex(x_val[example,...]))

# Plot
plt.figure(figsize=(11, 3), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(131)
plt.imshow(np.abs(predicted_image),cmap='gray',vmin=0,vmax=1)
plt.axis('off')
plt.colorbar()
plt.title('Predicted')

plt.subplot(132)
plt.imshow(np.abs(act_image),cmap='gray',vmin=0,vmax=1)
plt.axis('off')
plt.colorbar()
plt.title('True Image')

plt.subplot(133)
plt.imshow(np.abs(act_image-predicted_image),cmap='gray',vmin=0)
plt.axis('off')
plt.colorbar()
plt.title('Difference Image')

plt.show()<jupyter_output><empty_output><jupyter_text>Compare to least squares solution with dataHere we compare to an alterantive aproach, regularied least squares. In this technique, we build an encoding matrix which simulates the data acquisition. Then we minimize:$\parallel Ex-d \parallel_2 +  \lambda \parallel x \parallel_2$Where $\lambda$ is a factor that regularizes the solution when its illposed ( see https://en.wikipedia.org/wiki/Tikhonov_regularization ). The solution to this set of equations is:$ \widetilde{x} = (E^hE + \lambda I)^{-1}E^hd$Where I is an identity matrix. Similar to the neural network this is an aproximate solution.<jupyter_code># Lets also solve this a different way using a matrix inverse
def DFT_matrix(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    omega = np.exp( 2 * math.pi * 1J / N )
    W = np.power( omega, i * j ) / N #math.sqrt(N)
    return W

E = DFT_matrix(N)
E = np.fft.fftshift(E,axes=(0,))
E = E[idx,:]

# Grab the data
D = np.matrix.getH(channels_to_complex(kspace_val[400,...]))

# Solve for psuedo inverse
Eh = np.matrix.getH(E)
EhE = np.matmul(Eh,E)
Ei = np.linalg.inv(EhE + 0.000001*np.identity(N))
EiEh = np.matmul(Ei,Eh)

linear_algebra_prediction = np.transpose(np.matmul(EiEh,D))
    
plt.figure(figsize=(11, 11), dpi=80, facecolor='w', edgecolor='k')

plt.subplot(231)
plt.imshow(np.abs(linear_algebra_prediction),cmap='gray',vmin=0)
plt.axis('off')
plt.title('Least Squares Solution')
plt.subplot(234)
plt.imshow(np.abs(linear_algebra_prediction-act_image),cmap='gray',vmin=0,vmax=0.2)
plt.axis('off')
plt.title('Difference Least Squares')

plt.subplot(232)
plt.imshow(np.abs(predicted_image),cmap='gray',vmin=0,vmax=1)
plt.axis('off')
plt.title('Neural Net Prediction')
plt.subplot(235)
plt.imshow(np.abs(predicted_image-act_image),cmap='gray',vmin=0,vmax=0.2)
plt.axis('off')
plt.title('Difference Neural Net')

plt.subplot(233)
plt.imshow(np.abs(act_image),cmap='gray',vmin=0,vmax=1)
plt.axis('off')
plt.title('Actual Image')

plt.show()

print('Image Domain Mean Squared Error NN = ' + str(np.sum(np.square(abs(np.squeeze(predicted_image) - act_image)))) )
print('Image Domain Mean Squared Error LS = ' + str(np.sum(np.square(abs(linear_algebra_prediction - act_image)))) )

# Lets also get the kspace error
kspace_NN = np.matmul(E,np.squeeze(predicted_image))
kspace_LA = np.matmul(E,linear_algebra_prediction)

# Difference 
diff_kspace_NN = kspace_NN - D
diff_kspace_LA = kspace_LA - D
print('Kspace Mean Squared Error NN = ' + str(np.sum(np.square(abs(diff_kspace_NN)))) )
print('Kspace Mean Squared Error LS = ' + str(np.sum(np.square(abs(diff_kspace_LA)))) )<jupyter_output><empty_output><jupyter_text>Load real MRI data to testThis is actual acquired MRI data from a brain scan consisting. The data size is larger and we crop in k-space. Just to make things doable in a short time we are keeping everything 1D, as above.<jupyter_code># Load a Kspace dataset from an actual acquisition
with h5py.File('/content/drive/MyDrive/ML4MI_BOOTCAMP_DATA/Example_MRI_Data.h5','r') as hf:
    kspace_mri = np.array(hf['Kspace'])

#Crop Kspace
crop = ( kspace_mri.shape[-2] - N ) // 2
kspace_mri = kspace_mri[...,::2,crop:-crop]
    
print(f'Kspace size = {kspace_mri.shape} [ channels, slices, Nx, Ny], type = {kspace_mri.dtype}')
coils = kspace_mri.shape[0]
slices = kspace_mri.shape[1]<jupyter_output><empty_output><jupyter_text>Run a traditional reconstruction The most common reconstruction on MRI scanners is to just do a discrete Fourier transform of the data. Just a note, the data actually has 48 recievers of the signal. We are taking the sum of squares to average these signals.<jupyter_code># Traditional recon of fully sampled data
image_full = np.fft.ifftn(kspace_mri,axes=(-1,))

# do sum of squares to average coils (detectors)
image_full = np.sum(abs(image_full),axis=0)
image_full = np.sqrt(image_full)

# Make a montage (there are other options)
plot_image = montage(image_full[8::2,:,:])  
    
# Show the image
plt.figure(figsize=(20,20))
plt.imshow(plot_image,aspect=1,interpolation='bilinear',cmap='gray')
plt.axis('off')
plt.title('DFT of Kspace')
plt.show()<jupyter_output><empty_output><jupyter_text>Do inference on the real MRI data Machine Learning Based Reconstruction<jupyter_code># Subsample kspace and convert to channels
kspace_mri2 = kspace_mri[:,:,:,sampling_mask]
kspace_mri2 = np.stack((kspace_mri2.real,kspace_mri2.imag),axis=-1)


kspace_mri2 = np.reshape(kspace_mri2,(-1,N,number_phase_encodes,2))
print(kspace_mri2.shape)

# Run model
kspace_torch = torch.tensor(np.moveaxis( kspace_mri2, -1, 1))
predicted_image = recon_model(kspace_torch.to(device))
predicted_image = predicted_image.detach().cpu().numpy()
image_NN = np.moveaxis(predicted_image,1,-1)
print(image_NN.shape)

# Reshape
image_NN = np.reshape( image_NN,(coils,slices,N,N,2))
image_NN = channels_to_complex(image_NN)

# do sum of squares to average coils (detectors)
image_NN = np.sum(abs(image_NN),axis=0)
image_NN = np.sqrt(image_NN)

# Make a montage (there are other options)
plot_image = montage( image_NN[8::2,:,:])

# Show the image
plt.figure(figsize=(20,20))
plt.imshow(plot_image,aspect=1,interpolation='bilinear',cmap='gray')
plt.axis('off')
plt.title('Neural network prediction from Kspace')
plt.show()<jupyter_output><empty_output><jupyter_text>Linear algebra based solution<jupyter_code>image_LA = np.zeros(image_full.shape,dtype=image_full.dtype)

for k in range(slices):

  # Subsample kspace and convert to channels
  kspace_mri2 = np.squeeze(kspace_mri[:,k,:,:])
  kspace_mri2 = kspace_mri2[:,:,sampling_mask]

  kspace_mri2 = np.reshape(kspace_mri2,(-1,number_phase_encodes))
  kspace_mri2 = np.expand_dims(kspace_mri2,-1)
  
  # Also do for Least squares estimate
  image = np.matmul(EiEh,kspace_mri2)
  image = np.reshape(image,newshape=(coils,N,N))

  # do sum of squares to average coils (detectors)
  image = np.sum(abs(image),axis=0)
  image = np.sqrt(image)

  image_LA[k,:,:] = np.fliplr(image)

# Make a montage (there are other options)
plot_image = montage( image_LA[8::2,:,:])

# Show the image
plt.figure(figsize=(20,20))
plt.imshow(plot_image,aspect=1,interpolation='bilinear',cmap='gray')
plt.axis('off')
plt.title('Linear algebra prediction from Kspace')
plt.show()<jupyter_output><empty_output><jupyter_text>Now compare the solutions<jupyter_code>slice = 24

print(image_LA.shape)
print(image_NN.shape)
print(image_full.shape)

plt.figure(figsize=(20,20))
plt.subplot(131)
plt.imshow(abs(image_LA[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Linear Algebra')

plt.subplot(132)
plt.imshow(abs(image_NN[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Neural Net')

plt.subplot(133)
plt.imshow(abs(image_full[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Ground Truth')

plt.show()
# Slice for viewing
slice = 24

# Scale to minimize difference (scaling unimportant in MRI)
scale_LA = np.sum( image_full*np.conj(image_LA)) /np.sum(image_LA**2)
scale_NN = np.sum( image_full*np.conj(image_NN)) /np.sum(image_NN**2)

diff_LA = scale_LA*image_LA - image_full
diff_NN = scale_NN*image_NN - image_full

# Print Error
error_LA = np.linalg.norm(diff_LA)/np.linalg.norm(image_full)
error_NN = np.linalg.norm(diff_NN)/np.linalg.norm(image_full)

print(f'Image MSE Linear Algebra = {error_LA}')
print(f'Image MSE Neural Network = {error_NN}')

plt.figure(figsize=(20,20))
plt.subplot(131)
plt.imshow(abs(diff_LA[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Linear Algebra')

plt.subplot(132)
plt.imshow(abs(diff_NN[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Neural Net')

plt.subplot(133)
plt.imshow(abs(image_full[slice,:,:]),cmap='gray')
plt.axis('off')
plt.title('Ground Truth')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11990844742084517, 'fraction_numerical': 0.015004450472597804, 'mean_word_length': 4.967121901871523, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 53, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '581089', 'license': 'MIT', 'n_tokens_mistral': 8570, 'n_tokens_neox': 8136, 'n_words': 2188, 'path': 'ImageReconstruction/CoLab_AutoMap_Recon_Torch.ipynb', 'repo_name': 'kmjohnson3/ML4MI_BootCamp', 'size': '23593'}