diff --git a/demo/.demo_3D_shepp_logan.py.swp b/demo/.demo_3D_shepp_logan.py.swp
new file mode 100644
index 0000000000000000000000000000000000000000..47d71ec459aaf1b5b32aa615cd794a6374c1b342
Binary files /dev/null and b/demo/.demo_3D_shepp_logan.py.swp differ
diff --git a/demo/demo_3D_shepp_logan.py b/demo/demo_3D_shepp_logan.py
index 920aae8a5c6dacb2c0b0edf47ab4437ade04dcd1..e0cf4009c5e98df7ac627969ea0b01ee40b7ae2e 100644
--- a/demo/demo_3D_shepp_logan.py
+++ b/demo/demo_3D_shepp_logan.py
@@ -4,46 +4,32 @@ import mbircone
from demo_utils import plot_image, plot_gif, plt_cmp_3dobj
# Set sinogram shape
-num_det_rows = 200
+num_det_rows = 128
num_det_channels = 128
num_views = 144
# Reconstruction parameters
-sharpness = 0.2
-snr_db = 31.0
-
-# magnification is unitless.
-magnification = 2
-
-# All distances are in unit of 1 ALU = 1 mm.
-dist_source_detector = 600
-delta_pixel_detector = 0.9
-delta_pixel_image = 1
-channel_offset = 0
-row_offset = 0
-max_iterations = 20
+sharpness = 1.0 # Controls sharpness of reconstruction
+magnification = 2.0 # Ratio of (source to detector)/(source to center of rotation)
+dist_source_detector = 2*num_det_channels # distance from source to detector in ALU
+angles = np.linspace(0, 2 * np.pi, num_views, endpoint=False) # set of projection angles
# Display parameters
-vmin = 1.0
-vmax = 1.1
+vmin = 1.0 # minimum display value
+vmax = 1.1 # maximum display value
filename = 'output/3D_shepp_logan/results.png'
# Generate a 3D shepp logan phantom.
-ROR, boundary_size= mbircone.cone3D.compute_img_size(num_views, num_det_rows, num_det_channels,
- dist_source_detector,
- magnification,
- channel_offset=channel_offset, row_offset=row_offset,
- delta_pixel_detector=delta_pixel_detector,
- delta_pixel_image=delta_pixel_image)
-Nz, Nx, Ny = ROR
+ROR, boundary_size = mbircone.cone3D.compute_img_size(num_views, num_det_rows, num_det_channels, dist_source_detector, magnification)
+Nz, Nx, Ny = ROR # function ensures that Nx=Ny
img_slices_boundary_size, img_rows_boundary_size, img_cols_boundary_size = boundary_size
print('ROR of the recon is:', (Nz, Nx, Ny))
-# Set phantom parameters to generate a phantom inside ROI according to ROR and boundary_size.
-# All valid pixels should be inside ROI.
-num_rows_cols = Nx - 2 * img_rows_boundary_size # Assumes a square image
-num_slices_phantom = Nz - 2 * img_slices_boundary_size
-print('ROI and shape of phantom is:', num_slices_phantom, num_rows_cols, num_rows_cols)
+# Determine required phantom size based on ROI.
+# ROI is determine by ROR - 2*boundary_size
+num_rows_cols = Nx - 2 * img_rows_boundary_size # determines the width and height of the ROI
+num_slices_phantom = Nz - 2 * img_slices_boundary_size # determines the depth of the ROI
+print('ROI and shape of phantom is:', num_slices_phantom, num_rows_cols, num_rows_cols) # assumes that width=height
# Set display indexes
display_slice = img_slices_boundary_size + int(0.4*num_slices_phantom)
@@ -51,22 +37,17 @@ display_x = num_rows_cols // 2
display_y = num_rows_cols // 2
display_view = 0
-# Generate a phantom
-phantom = mbircone.phantom.gen_shepp_logan_3d(num_rows_cols, num_rows_cols, num_slices_phantom)
-print('Generated phantom shape = ', np.shape(phantom))
-phantom = mbircone.cone3D.pad_roi2ror(phantom, boundary_size)
-print('Padded phantom shape = ', np.shape(phantom))
+# Generate a 3D shepp logan phantom
+phantom = mbircone.phantom.gen_shepp_logan_3d(num_rows_cols, num_rows_cols, num_slices_phantom) # generate phatom within ROI
+print('Generated ROI phantom shape = ', np.shape(phantom))
+phantom = mbircone.cone3D.pad_roi2ror(phantom, boundary_size) # generate phantom wihtin ROR
+print('Padded ROR phantom shape = ', np.shape(phantom))
-# Generate simulated data using forward projector on the 3D shepp logan phantom.
-angles = np.linspace(0, 2 * np.pi, num_views, endpoint=False)
+# Compute the forward projection of the phantom
+sino = mbircone.cone3D.project(phantom, angles, num_det_rows, num_det_channels, dist_source_detector, magnification)
-# After setting the geometric parameter, the shape of the input phantom should be equal to the calculated geometric parameter.
-# Input a phantom with wrong shape will generate a bunch of issue in C.
-sino = mbircone.cone3D.project(phantom, angles,
- num_det_rows=num_det_rows, num_det_channels=num_det_channels,
- dist_source_detector=dist_source_detector, magnification=magnification,
- delta_pixel_detector=delta_pixel_detector, delta_pixel_image=delta_pixel_image,
- channel_offset=channel_offset, row_offset=row_offset)
+# Reconstruction sinogram data
+recon = mbircone.cone3D.recon(sino, angles, dist_source_detector, magnification, sharpness=sharpness )
# create output folder
os.makedirs('output/3D_shepp_logan/', exist_ok=True)
@@ -75,12 +56,6 @@ print('sino shape = ', np.shape(sino), sino.dtype)
plot_image(sino[display_view, :, :], title='sino', filename='output/3D_shepp_logan/sino-shepp-logan-3D-view(%.2f).png' % angles[0])
# plot_gif(sino, 'output', 'sino-shepp-logan-3D')
-recon = mbircone.cone3D.recon(sino, angles,
- dist_source_detector=dist_source_detector, magnification=magnification,
- delta_pixel_detector=delta_pixel_detector, delta_pixel_image=delta_pixel_image,
- channel_offset=channel_offset, row_offset=row_offset,
- sharpness=sharpness, snr_db=snr_db, max_iterations=max_iterations)
-
print('recon shape = ', np.shape(recon))
np.save('output/3D_shepp_logan/recon.npy', recon)
diff --git a/demo/demo_mace3D.py b/demo/demo_mace3D.py
index 4e4bf05e6c5e00581ab9e802aef5bea085aae16b..f0cdc371bae99ccd373f03d1b34a4d24e6a64042 100644
--- a/demo/demo_mace3D.py
+++ b/demo/demo_mace3D.py
@@ -3,22 +3,21 @@ import numpy as np
import math
import urllib.request
import tarfile
-from keras.models import model_from_json
import mbircone
import demo_utils, denoiser_utils
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
"""
-This script is a demonstration of the mace3D reconstruction algorithm. Demo functionality includes
+This script is a demonstration of the mace3D reconstruction algorithm. Demo functionality includes
* downloading phantom and denoiser data from specified urls
- * generating sinogram by projecting the phantom and then adding transmission noise
- * performing a 3D MACE reconstruction.
+ * generating synthetic sinogram data by forward projecting the phantom and then adding transmission noise
+ * performing a 3D MACE and qGGMRF reconstructions and displaying them.
"""
-print('This script is a demonstration of the mace3D reconstruction algorithm. Demo functionality includes \
+print('This script is a demonstration of the mace3D reconstruction algorithm. Demo functionality includes \
\n\t * downloading phantom and denoiser data from specified urls \
-\n\t * generating sinogram by projecting the phantom and then adding transmission noise\
-\n\t * performing a 3D MACE reconstruction.')
+\n\t * generating synthetic sinogram data by forward projecting the phantom and then adding transmission noise\
+\n\t * performing a 3D MACE and qGGMRF reconstructions and displaying them.')
# ###########################################################################
# Set the parameters to get the data and do the recon
@@ -26,11 +25,6 @@ print('This script is a demonstration of the mace3D reconstruction algorithm. D
# Change the parameters below for your own use case.
-# Denoiser function to be used in MACE. For the built-in demo, this should be one of dncnn_keras or dncnn_ct
-# Other denoisers built in keras can be used with minimal modification by putting the architecture and weights
-# in model.json and model.hdf5 in the denoiser_path set below
-denoiser_type = 'dncnn_ct'
-
# The url to the data repo.
data_repo_url = 'https://github.com/cabouman/mbir_data/raw/master/'
@@ -38,11 +32,10 @@ data_repo_url = 'https://github.com/cabouman/mbir_data/raw/master/'
# This file will be used to retrieve urls to files that we are going to download
yaml_url = os.path.join(data_repo_url, 'index.yaml')
-# Choice of phantom and denoiser files.
+# Choice of phantom file.
# These should be valid choices specified in the index file.
-# The urls to phantom data and NN weights will be parsed from data_repo_url and the choices of phantom and denoiser specified below.
+# The url to phantom data will be parsed from data_repo_url and the choices of phantom specified below.
phantom_name = 'bottle_cap_3D'
-denoiser_name = denoiser_type
# destination path to download and extract the phantom and NN weight files.
target_dir = './demo_data/'
@@ -62,9 +55,8 @@ num_views = 75 # number of projection views
sino_noise_sigma = 0.01 # transmission noise level
# MACE recon parameters
+sharpness = 1.0 # Parameter to control regularization level of reconstruction.
max_admm_itr = 10 # max ADMM iterations for MACE reconstruction
-sharpness = 1.0
-prior_weight = 0.5 # cumulative weights for three prior agents.
# ######### End of parameters #########
@@ -78,19 +70,13 @@ index_path = demo_utils.download_and_extract(yaml_url, target_dir)
url_index = demo_utils.load_yaml(index_path)
# get urls to phantom and denoiser parameter file
phantom_url = os.path.join(data_repo_url, url_index['phantom'][phantom_name]) # url to download the 3D image volume phantom file
-denoiser_url = os.path.join(data_repo_url, url_index['denoiser'][denoiser_name]) # url to download the denoiser parameter file
+denoiser_url = os.path.join(data_repo_url, url_index['denoiser']['dncnn_ct']) # url to download the denoiser parameter file
# download phantom file
phantom_path = demo_utils.download_and_extract(phantom_url, target_dir)
-# download and extract NN weights and structure files
+# download and extract NN weights and structure files used for MACE denoiser
denoiser_path = demo_utils.download_and_extract(denoiser_url, target_dir)
-
-# ###########################################################################
-# Generate downsampled phantom
-# ###########################################################################
-print("Generating downsampled 3D phantom volume ...")
-
# load original phantom
phantom = np.load(phantom_path)
print("shape of phantom = ", phantom.shape)
@@ -107,54 +93,38 @@ sino = mbircone.cone3D.project(phantom, angles,
num_det_rows, num_det_channels,
dist_source_detector, magnification,
delta_pixel_detector=delta_pixel_detector)
-weights = mbircone.cone3D.calc_weights(sino, weight_type='transmission')
+sino_weights = mbircone.cone3D.calc_weights(sino, weight_type='transmission')
# Add transmission noise
-noise = sino_noise_sigma * 1. / np.sqrt(weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))
+noise = sino_noise_sigma * 1. / np.sqrt(sino_weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))
sino_noisy = sino + noise
+# ###########################################################################
+# Perform qGGMRF reconstruction
+# ###########################################################################
+print("Performing qGGMRF reconstruction ...")
+recon_qGGMRF = mbircone.cone3D.recon(sino_noisy, angles, dist_source_detector, magnification,
+ delta_pixel_detector=delta_pixel_detector,
+ sharpness=sharpness, weight_type='transmission',
+ verbose=1)
# ###########################################################################
# Set up the denoiser
# ###########################################################################
-# This demo includes a custom CNN trained on CT images and a generic DnCNN in keras
-# The choice is set in the parameter section above.
+# The MACE reconstruction algorithm requires the use of a denoiser as a prior model.
+# This demo includes a custom DnCNN denoiser trained on CT images, which can be used for other applications.
print("Loading denoiser function and model ...")
-
-# DnCNN denoiser trained on CT images. This is the denoiser that we recommend using.
-if denoiser_type == 'dncnn_ct':
- print("Denoiser function: custom DnCNN trained on CT images.")
-
- # Load denoiser model structure and weights
- denoiser_model = denoiser_utils.DenoiserCT(checkpoint_dir=os.path.join(denoiser_path, 'model_dncnn_ct'))
- # Define the denoiser using this model. This version requires some interface code to match with MACE.
- def denoiser(img_noisy):
- img_noisy = np.expand_dims(img_noisy, axis=1)
- upper_range = denoiser_utils.calc_upper_range(img_noisy)
- img_noisy = img_noisy/upper_range
- testData_obj = denoiser_utils.DataLoader(img_noisy)
- img_denoised = denoiser_model.denoise(testData_obj, batch_size=256)
- img_denoised = img_denoised*upper_range
- return np.squeeze(img_denoised)
-# DnCNN denoiser in Keras. This denoiser model is trained on natural images.
-elif denoiser_type == 'dncnn_keras':
- print("Denoiser function: DnCNN trained on natural images.")
-
- # Load denoiser model structure and weights
- json_path = os.path.join(denoiser_path, 'model_dncnn_keras/model.json') # model architecture file
- weight_path = os.path.join(denoiser_path, 'model_dncnn_keras/model.hdf5') # model weight file
- with open(json_path, 'r') as json_file:
- denoiser_model = model_from_json(json_file.read()) # load model architecture
-
- denoiser_model.load_weights(weight_path) # load model weights
-
- # Define the denoiser using this model.
- def denoiser(img_noisy):
- img_noisy = img_noisy[..., np.newaxis] # (Nz,N0,N1,...,Nm,1)
- img_denoised = denoiser_model.predict(img_noisy) # inference
- return np.squeeze(img_denoised)
-else:
- raise RuntimeError('Unkown denoiser_type. Should be either dncnn_ct or dncnn_keras.')
+# Load denoiser model structure and pre-trained model weights
+denoiser_model = denoiser_utils.DenoiserCT(checkpoint_dir=os.path.join(denoiser_path, 'model_dncnn_ct'))
+# Define the denoiser using this model. This version requires some interface code to match with MACE.
+def denoiser(img_noisy):
+ img_noisy = np.expand_dims(img_noisy, axis=1)
+ upper_range = denoiser_utils.calc_upper_range(img_noisy)
+ img_noisy = img_noisy/upper_range
+ testData_obj = denoiser_utils.DataLoader(img_noisy)
+ img_denoised = denoiser_model.denoise(testData_obj, batch_size=256)
+ img_denoised = img_denoised*upper_range
+ return np.squeeze(img_denoised)
# ###########################################################################
@@ -163,26 +133,22 @@ else:
print("Performing MACE reconstruction ...")
recon_mace = mbircone.mace.mace3D(sino_noisy, angles, dist_source_detector, magnification,
denoiser=denoiser, denoiser_args=(),
- max_admm_itr=max_admm_itr, prior_weight=prior_weight,
+ max_admm_itr=max_admm_itr,
+ init_image=recon_qGGMRF,
+ sharpness=sharpness,
delta_pixel_detector=delta_pixel_detector,
weight_type='transmission',
- sharpness=sharpness,
- save_path=save_path)
+ verbose=1)
recon_shape = recon_mace.shape
print("Reconstruction shape = ", recon_shape)
# ###########################################################################
-# Post-process reconstruction results
+# Generating phantom and reconstruction images
# ###########################################################################
-print("Post processing MACE reconstruction results ...")
-# Save recon results as a numpy array
-np.save(os.path.join(save_path, "recon_mace.npy"), recon_mace)
-# load qGGMRF recon
-recon_qGGMRF = np.load(os.path.join(save_path, "recon_qGGMRF.npy"))
-
+print("Generating phantom and reconstruction images ...")
# Plot axial slices of phantom and recon
-display_slices = [7, 10, 13, 16, 19, 22]
+display_slices = [7, 12, 17, 22]
for display_slice in display_slices:
demo_utils.plot_image(phantom[display_slice], title=f'phantom, axial slice {display_slice}',
filename=os.path.join(save_path, f'phantom_slice{display_slice}.png'), vmin=0, vmax=0.5)
@@ -190,7 +156,6 @@ for display_slice in display_slices:
filename=os.path.join(save_path, f'recon_mace_slice{display_slice}.png'), vmin=0, vmax=0.5)
demo_utils.plot_image(recon_qGGMRF[display_slice], title=f'qGGMRF reconstruction, axial slice {display_slice}',
filename=os.path.join(save_path, f'recon_qGGMRF_slice{display_slice}.png'), vmin=0, vmax=0.5)
-
# Plot 3D phantom and recon image volumes as gif images.
demo_utils.plot_gif(phantom, save_path, 'phantom', vmin=0, vmax=0.5)
demo_utils.plot_gif(recon_mace, save_path, 'recon_mace', vmin=0, vmax=0.5)
diff --git a/demo/demo_mace3D_fast.py b/demo/demo_mace3D_fast.py
index dbdc46dddb1f3777912ea82501405799a22c37e2..7b7a8f86ae48a2530a3a88de0ac98ffddc474b95 100644
--- a/demo/demo_mace3D_fast.py
+++ b/demo/demo_mace3D_fast.py
@@ -3,24 +3,21 @@ import numpy as np
import math
import urllib.request
import tarfile
-from keras.models import model_from_json
import mbircone
import demo_utils, denoiser_utils
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
"""
-This script is a quick demonstration of the mace3D reconstruction algorithm. Demo functionality includes
+This script is a fast demonstration of the mace3D reconstruction algorithm that uses a low-res phantom. Demo functionality includes
* downloading phantom and denoiser data from specified urls
- * downsampling the phantom along all three dimensions
- * generating sinogram by projecting the phantom and then adding transmission noise
- * performing a 3D MACE reconstruction.
+ * generating synthetic sinogram data by forward projecting the phantom and then adding transmission noise
+ * performing a 3D MACE and qGGMRF reconstructions and displaying them.
"""
-print('This script is a quick demonstration of the mace3D reconstruction algorithm. Demo functionality includes \
+print('This script is a demonstration of the mace3D reconstruction algorithm. Demo functionality includes \
\n\t * downloading phantom and denoiser data from specified urls \
-\n\t * downsampling the phantom along all three dimensions \
-\n\t * generating sinogram by projecting the phantom and then adding transmission noise\
-\n\t * performing a 3D MACE reconstruction.')
+\n\t * generating synthetic sinogram data by forward projecting the phantom and then adding transmission noise\
+\n\t * performing a 3D MACE and qGGMRF reconstructions and displaying them.')
# ###########################################################################
# Set the parameters to get the data and do the recon
@@ -28,11 +25,6 @@ print('This script is a quick demonstration of the mace3D reconstruction algorit
# Change the parameters below for your own use case.
-# Denoiser function to be used in MACE. For the built-in demo, this should be one of dncnn_keras or dncnn_ct
-# Other denoisers built in keras can be used with minimal modification by putting the architecture and weights
-# in model.json and model.hdf5 in the denoiser_path set below
-denoiser_type = 'dncnn_ct'
-
# The url to the data repo.
data_repo_url = 'https://github.com/cabouman/mbir_data/raw/master/'
@@ -40,15 +32,14 @@ data_repo_url = 'https://github.com/cabouman/mbir_data/raw/master/'
# This file will be used to retrieve urls to files that we are going to download
yaml_url = os.path.join(data_repo_url, 'index.yaml')
-# Choice of phantom and denoiser files.
+# Choice of phantom file.
# These should be valid choices specified in the index file.
-# The urls to phantom data and NN weights will be parsed from data_repo_url and the choices of phantom and denoiser specified below.
+# The url to phantom data will be parsed from data_repo_url and the choices of phantom specified below.
phantom_name = 'bottle_cap_3D'
-denoiser_name = denoiser_type
-# destination path to download and extract the phantom and NN weight files.
+# Destination path to download and extract the phantom and NN weight files.
target_dir = './demo_data/'
-# path to store output recon images
+# Path to store output recon images
save_path = './output/mace3D_fast/'
os.makedirs(save_path, exist_ok=True)
@@ -63,10 +54,9 @@ num_det_channels = 120 # number of detector channels
num_views = 50 # number of projection views
sino_noise_sigma = 0.005 # transmission noise level
-# MACE recon parameters
+# Recon parameters
+sharpness = 1.0 # Parameter to control regularization level of reconstruction.
max_admm_itr = 10 # max ADMM iterations for MACE reconstruction
-sharpness = 1.0
-prior_weight = 0.5 # cumulative weights for three prior agents.
# ######### End of parameters #########
@@ -80,7 +70,7 @@ index_path = demo_utils.download_and_extract(yaml_url, target_dir)
url_index = demo_utils.load_yaml(index_path)
# get urls to phantom and denoiser parameter file
phantom_url = os.path.join(data_repo_url, url_index['phantom'][phantom_name]) # url to download the 3D image volume phantom file
-denoiser_url = os.path.join(data_repo_url, url_index['denoiser'][denoiser_name]) # url to download the denoiser parameter file
+denoiser_url = os.path.join(data_repo_url, url_index['denoiser']['dncnn_ct']) # url to download the denoiser parameter file
# download phantom file
phantom_path = demo_utils.download_and_extract(phantom_url, target_dir)
@@ -116,56 +106,37 @@ sino = mbircone.cone3D.project(phantom, angles,
num_det_rows, num_det_channels,
dist_source_detector, magnification,
delta_pixel_detector=delta_pixel_detector)
-weights = mbircone.cone3D.calc_weights(sino, weight_type='transmission')
+sino_weights = mbircone.cone3D.calc_weights(sino, weight_type='transmission')
# Add transmission noise
-noise = sino_noise_sigma * 1. / np.sqrt(weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))
+noise = sino_noise_sigma * 1. / np.sqrt(sino_weights) * np.random.normal(size=(num_views, num_det_rows, num_det_channels))
sino_noisy = sino + noise
+# ###########################################################################
+# Perform qGGMRF reconstruction
+# ###########################################################################
+print("Performing qGGMRF reconstruction ...")
+recon_qGGMRF = mbircone.cone3D.recon(sino_noisy, angles, dist_source_detector, magnification,
+ delta_pixel_detector=delta_pixel_detector,
+ sharpness=sharpness, weight_type='transmission',
+ verbose=1)
# ###########################################################################
-# Set up the denoiser
+# Set up the denoiser used for mace3D recon
# ###########################################################################
-# This demo includes a custom CNN trained on CT images and a generic DnCNN in keras
-# The choice is set in the parameter section above.
+# This demo includes a custom DnCNN denoiser trained on CT images.
print("Loading denoiser function and model ...")
-
-# DnCNN denoiser trained on CT images. This is the denoiser that we recommend using.
-if denoiser_type == 'dncnn_ct':
- print("Denoiser function: custom DnCNN trained on CT images.")
-
- # Load denoiser model structure and weights
- denoiser_model = denoiser_utils.DenoiserCT(checkpoint_dir=os.path.join(denoiser_path, 'model_dncnn_ct'))
-
- # Define the denoiser using this model. This version requires some interface code to match with MACE.
- def denoiser(img_noisy):
- img_noisy = np.expand_dims(img_noisy, axis=1)
- upper_range = denoiser_utils.calc_upper_range(img_noisy)
- img_noisy = img_noisy/upper_range
- testData_obj = denoiser_utils.DataLoader(img_noisy)
- img_denoised = denoiser_model.denoise(testData_obj, batch_size=256)
- img_denoised = img_denoised*upper_range
- return np.squeeze(img_denoised)
-# DnCNN denoiser in Keras. This denoiser model is trained on natural images.
-elif denoiser_type == 'dncnn_keras':
- print("Denoiser function: DnCNN trained on natural images.")
-
- # Load denoiser model structure and weights
- json_path = os.path.join(denoiser_path, 'model_dncnn_keras/model.json') # model architecture file
- weight_path = os.path.join(denoiser_path, 'model_dncnn_keras/model.hdf5') # model weight file
- with open(json_path, 'r') as json_file:
- denoiser_model = model_from_json(json_file.read()) # load model architecture
-
- denoiser_model.load_weights(weight_path) # load model weights
-
- # Define the denoiser using this model.
- def denoiser(img_noisy):
- img_noisy = img_noisy[..., np.newaxis] # (Nz,N0,N1,...,Nm,1)
- img_denoised = denoiser_model.predict(img_noisy) # inference
- return np.squeeze(img_denoised)
-else:
- raise RuntimeError('Unkown denoiser_type. Should be either dncnn_ct or dncnn_keras.')
-
+# Load denoiser model structure and pre-trained model weights
+denoiser_model = denoiser_utils.DenoiserCT(checkpoint_dir=os.path.join(denoiser_path, 'model_dncnn_ct'))
+# Define the denoiser using this model. This version requires some interface code to match with MACE.
+def denoiser(img_noisy):
+ img_noisy = np.expand_dims(img_noisy, axis=1)
+ upper_range = denoiser_utils.calc_upper_range(img_noisy)
+ img_noisy = img_noisy/upper_range
+ testData_obj = denoiser_utils.DataLoader(img_noisy)
+ img_denoised = denoiser_model.denoise(testData_obj, batch_size=256)
+ img_denoised = img_denoised*upper_range
+ return np.squeeze(img_denoised)
# ###########################################################################
# Perform MACE reconstruction
@@ -173,26 +144,22 @@ else:
print("Performing MACE reconstruction ...")
recon_mace = mbircone.mace.mace3D(sino_noisy, angles, dist_source_detector, magnification,
denoiser=denoiser, denoiser_args=(),
- max_admm_itr=max_admm_itr, prior_weight=prior_weight,
+ max_admm_itr=max_admm_itr,
+ init_image=recon_qGGMRF,
sharpness=sharpness,
delta_pixel_detector=delta_pixel_detector,
weight_type='transmission',
- save_path=save_path)
+ verbose=1)
recon_shape = recon_mace.shape
print("Reconstruction shape = ", recon_shape)
# ###########################################################################
-# Post-process reconstruction results
+# Generating phantom and reconstruction images
# ###########################################################################
-print("Post processing MACE reconstruction results ...")
-# Save recon results as a numpy array
-np.save(os.path.join(save_path, "recon_mace.npy"), recon_mace)
-# load qGGMRF recon
-recon_qGGMRF = np.load(os.path.join(save_path, "recon_qGGMRF.npy"))
-
+print("Generating phantom and reconstruction images ...")
# Plot axial slices of phantom and recon
-display_slices = [7, 10, 13, 16, 19, 22]
+display_slices = [7, 12, 17, 22]
for display_slice in display_slices:
demo_utils.plot_image(phantom[display_slice], title=f'phantom, axial slice {display_slice}',
filename=os.path.join(save_path, f'phantom_slice{display_slice}.png'), vmin=0, vmax=0.5)
@@ -201,7 +168,7 @@ for display_slice in display_slices:
demo_utils.plot_image(recon_qGGMRF[display_slice], title=f'qGGMRF reconstruction, axial slice {display_slice}',
filename=os.path.join(save_path, f'recon_qGGMRF_slice{display_slice}.png'), vmin=0, vmax=0.5)
# Plot 3D phantom and recon image volumes as gif images.
-demo_utils.plot_gif(phantom, save_path, 'phantom_resized', vmin=0, vmax=0.5)
+demo_utils.plot_gif(phantom, save_path, 'phantom', vmin=0, vmax=0.5)
demo_utils.plot_gif(recon_mace, save_path, 'recon_mace', vmin=0, vmax=0.5)
demo_utils.plot_gif(recon_qGGMRF, save_path, 'recon_qGGMRF', vmin=0, vmax=0.5)
diff --git a/demo/demo_utils.py b/demo/demo_utils.py
index abf271d8c27bc3ee0f7654906da888ad3c183ff7..925bfadbe7155a139f87da56b1dbf81b084c83df 100644
--- a/demo/demo_utils.py
+++ b/demo/demo_utils.py
@@ -74,7 +74,7 @@ def plt_cmp_3dobj(phantom, recon, display_slice=None, display_x=None, display_y=
font_setting()
# display phantom
- fig, axs = plt.subplots(2, 3)
+ fig, axs = plt.subplots(2, 3, figsize=(20,10))
title = f'Phantom: Axial Scan {display_slice:d}.'
axs[0, 0].imshow(phantom[display_slice], vmin=vmin, vmax=vmax, cmap='gray', interpolation='none')
@@ -231,7 +231,7 @@ def download_and_extract(download_url, save_dir):
if save_path.endswith(('.tar', '.tar.gz')):
if is_download:
tar_file = tarfile.open(save_path)
- print("Extracting tarball file to {save_dir} ...")
+ print(f"Extracting tarball file to {save_dir} ...")
# Extract to save_dir.
tar_file.extractall(save_dir)
tar_file.close
diff --git a/mbircone/mace.py b/mbircone/mace.py
index a7acb171ad6928251eb9332eb3efe22f2bf3b535..d04fdbff8b25ada8aff0766448e40180d7aef804 100644
--- a/mbircone/mace.py
+++ b/mbircone/mace.py
@@ -53,16 +53,16 @@ def denoiser_wrapper(image_noisy, denoiser, denoiser_args, permute_vector, posit
def mace3D(sino, angles, dist_source_detector, magnification,
- denoiser, denoiser_args=(),
- max_admm_itr=10, rho=0.5, prior_weight=0.5,
- init_image=None, save_path=None,
- channel_offset=0.0, row_offset=0.0, rotation_offset=0.0,
- delta_pixel_detector=1.0, delta_pixel_image=None, ror_radius=None,
- sigma_y=None, snr_db=30.0, weights=None, weight_type='unweighted',
- positivity=True, p=1.2, q=2.0, T=1.0, num_neighbors=6,
- sharpness=0.0, sigma_x=None, sigma_p=None, max_iterations=3, stop_threshold=0.02,
- num_threads=None, NHICD=False, verbose=1,
- lib_path=__lib_path):
+ denoiser, denoiser_args=(),
+ max_admm_itr=10, rho=0.5, prior_weight=0.5,
+ init_image=None,
+ channel_offset=0.0, row_offset=0.0, rotation_offset=0.0,
+ delta_pixel_detector=1.0, delta_pixel_image=None, ror_radius=None,
+ sigma_y=None, snr_db=30.0, weights=None, weight_type='unweighted',
+ positivity=True, p=1.2, q=2.0, T=1.0, num_neighbors=6,
+ sharpness=0.0, sigma_x=None, sigma_p=None, max_iterations=3, stop_threshold=0.02,
+ num_threads=None, NHICD=False, verbose=1,
+ lib_path=__lib_path):
"""Computes 3-D conebeam beam reconstruction with multi-slice MACE alogorithm by fusing forward model proximal map with 2D denoisers across xy, xz, and yz planes.
Required arguments:
@@ -81,7 +81,6 @@ def mace3D(sino, angles, dist_source_detector, magnification,
- **rho** (*float*): [Default=0.5] step size of ADMM update in MACE, range (0,1). The value of ``rho`` mainly controls the convergence speed of MACE algorithm.
- **prior_weight** (*ndarray*): [Default=0.5] weights for prior agents, specified by either a scalar value or a 1D array. If a scalar is specified, then all three prior agents use the same weight of (prior_weight/3). If an array is provided, then the array should have three elements corresponding to the weight of denoisers in XY, YZ, and XZ planes respectively. The weight for forward model proximal map agent will be calculated as 1-sum(prior_weight) so that the sum of all agent weights are equal to 1. Each entry of prior_weight should have value between 0 and 1. sum(prior_weight) needs to be no greater than 1.
- **init_image** (*ndarray, optional*): [Default=None] Initial value of MACE reconstruction image, specified by either a scalar value or a 3-D numpy array with shape (num_img_slices,num_img_rows,num_img_cols). If None, the inital value of MACE will be automatically determined by a qGGMRF reconstruction.
- - **save_path** (*str, optional*): [Default=None] Path to directory that saves the intermediate results of MACE.
Optional arguments inherited from ``cone3D.recon`` (with changing default value of ``max_iterations``):
- **channel_offset** (*float, optional*): [Default=0.0] Distance in :math:`ALU` from center of detector to the source-detector line along a row.
- **row_offset** (*float, optional*): [Default=0.0] Distance in :math:`ALU` from center of detector to the source-detector line along a column.
@@ -133,10 +132,14 @@ def mace3D(sino, angles, dist_source_detector, magnification,
# Calculate automatic value of sigma_p
if sigma_p is None:
sigma_p = cone3D.auto_sigma_p(sino, delta_pixel_detector=delta_pixel_detector, sharpness=sharpness)
+ # if no init_image is provided, then use qGGMRF recon as init_image.
if init_image is None:
if verbose:
start = time.time()
print("Computing qGGMRF reconstruction. This will be used as MACE initialization point.")
+ if sigma_x is None:
+ sigma_x = cone3D.auto_sigma_x(sino, delta_pixel_detector=delta_pixel_detector, sharpness=sharpness)
+ # qGGMRF recon
init_image = cone3D.recon(sino, angles, dist_source_detector, magnification,
channel_offset=channel_offset, row_offset=row_offset, rotation_offset=rotation_offset,
delta_pixel_detector=delta_pixel_detector, delta_pixel_image=delta_pixel_image, ror_radius=ror_radius,
@@ -148,9 +151,6 @@ def mace3D(sino, angles, dist_source_detector, magnification,
end = time.time()
elapsed_t = end-start
print(f"Done computing qGGMRF reconstruction. Elapsed time: {elapsed_t:.2f} sec.")
- if not (save_path is None):
- print("Save qGGMRF reconstruction to disk.")
- np.save(os.path.join(save_path, 'recon_qGGMRF.npy'), init_image)
if np.isscalar(init_image):
(num_views, num_det_rows, num_det_channels) = np.shape(sino)
@@ -211,11 +211,6 @@ def mace3D(sino, angles, dist_source_detector, magnification,
if verbose:
denoise_elapsed = denoise_end - denoise_start
print(f"Done denoising in all hyper-planes, elapsed time {denoise_elapsed:.2f} sec")
- if not (save_path is None):
- for i in range(4):
- np.save(os.path.join(save_path, f'X{i}_itr{itr}.npy'), X[i])
- np.save(os.path.join(save_path, f'W{i}_itr{itr}.npy'), W[i])
-
Z = sum([beta[k]*(2*X[k]-W[k]) for k in range(image_dim+1)])
for k in range(image_dim+1):
W[k] += 2*rho*(Z-X[k])
@@ -232,7 +227,7 @@ def mace3D(sino, angles, dist_source_detector, magnification,
def mace4D(sino, angles, dist_source_detector, magnification,
denoiser, denoiser_args=(),
max_admm_itr=10, rho=0.5, prior_weight=0.5,
- init_image=None, save_path=None,
+ init_image=None,
cluster_ticket=None,
channel_offset=0.0, row_offset=0.0, rotation_offset=0.0,
delta_pixel_detector=1.0, delta_pixel_image=None, ror_radius=None,
@@ -266,7 +261,6 @@ def mace4D(sino, angles, dist_source_detector, magnification,
- **rho** (*float*): [Default=0.5] step size of ADMM update in MACE, range (0,1). The value of ``rho`` mainly controls the convergence speed of MACE algorithm.
- **prior_weight** (*ndarray*): [Default=0.5] weights for prior agents, specified by either a scalar value or a 1D array. If a scalar is specified, then all three prior agents use the same weight of (prior_weight/3). If an array is provided, then the array should have three elements corresponding to the weight of denoisers in XY-t, YZ-t, and XZ-t hyperplanes respectively. The weight for forward model proximal map agent will be calculated as 1-sum(prior_weight) so that the sum of all agent weights are equal to 1. Each entry of prior_weight should have value between 0 and 1. sum(prior_weight) needs to be no greater than 1.
- **init_image** (*ndarray, optional*): [Default=None] Initial value of MACE reconstruction image, specified by either a scalar value, a 3-D numpy array with shape (num_img_slices,num_img_rows,num_img_cols), or a 4-D numpy array with shape (num_time_points, num_img_slices,num_img_rows,num_img_cols). If None, the inital value of MACE will be automatically determined by a stack of 3-D qGGMRF reconstructions at different time points.
- - **save_path** (*str, optional*): [Default=None] Path to directory that saves the intermediate results of MACE. If not None, the inital qGGMRF reconstruction and input/output images of all agents from each MACE iteration will be saved to ``save_path``.
Arguments specific to multi-node computation:
- **cluster_ticket** (*Object*): [Default=None] A ticket used to access a specific cluster, that can be obtained from ``multinode.get_cluster_ticket``. If cluster_ticket is None, the process will run in serial. See `dask_jobqueue <https://jobqueue.dask.org/en/latest/api.html>`_ for more information.
Optional arguments inherited from ``cone3D.recon`` (with changing default value of ``max_iterations``):
@@ -340,25 +334,29 @@ def mace4D(sino, angles, dist_source_detector, magnification,
# Calculate automatic value of sigma_p
if sigma_p is None:
sigma_p = cone3D.auto_sigma_p(sino, delta_pixel_detector=delta_pixel_detector, sharpness=sharpness)
- # Fixed args dictionary used for multi-node parallelization
- constant_args = {'dist_source_detector':dist_source_detector, 'magnification':magnification,
- 'channel_offset':channel_offset, 'row_offset':row_offset, 'rotation_offset':rotation_offset,
- 'delta_pixel_detector':delta_pixel_detector, 'delta_pixel_image':delta_pixel_image, 'ror_radius':ror_radius,
- 'sigma_y':sigma_y, 'sigma_p':sigma_p, 'sigma_x':sigma_x,
- 'positivity':positivity, 'p':p, 'q':q, 'T':T, 'num_neighbors':num_neighbors,
- 'max_iterations':20, 'stop_threshold':stop_threshold,
- 'num_threads':num_threads, 'NHICD':NHICD, 'verbose':qGGMRF_verbose, 'lib_path':lib_path
- }
- # List of variable args dictionaries used for multi-node parallelization
- variable_args_list = [{'sino': sino[t], 'angles':angles[t], 'weights':weights[t]}
- for t in range(Nt)]
+ if cluster_ticket is not None:
+ # Fixed args dictionary used for multi-node parallelization
+ constant_args = {'dist_source_detector':dist_source_detector, 'magnification':magnification,
+ 'channel_offset':channel_offset, 'row_offset':row_offset, 'rotation_offset':rotation_offset,
+ 'delta_pixel_detector':delta_pixel_detector, 'delta_pixel_image':delta_pixel_image, 'ror_radius':ror_radius,
+ 'sigma_y':sigma_y, 'sigma_p':sigma_p,
+ 'positivity':positivity, 'p':p, 'q':q, 'T':T, 'num_neighbors':num_neighbors,
+ 'max_iterations':20, 'stop_threshold':stop_threshold,
+ 'num_threads':num_threads, 'NHICD':NHICD, 'verbose':qGGMRF_verbose, 'lib_path':lib_path
+ }
+ # List of variable args dictionaries used for multi-node parallelization
+ variable_args_list = [{'sino': sino[t], 'angles':angles[t], 'weights':weights[t]}
+ for t in range(Nt)]
# if init_image is not provided, use qGGMRF recon result as init_image.
if init_image is None:
if verbose:
start = time.time()
print("Computing qGGMRF reconstruction at all time points. This will be used as MACE initialization point.")
-
- if not (cluster_ticket is None):
+ if sigma_x is None:
+ sigma_x = cone3D.auto_sigma_x(sino, delta_pixel_detector=delta_pixel_detector, sharpness=sharpness)
+
+ if cluster_ticket is not None:
+ constant_args['sigma_x'] = sigma_x
init_image = np.array(multinode.scatter_gather(cluster_ticket, cone3D.recon,
variable_args_list=variable_args_list,
constant_args=constant_args,
@@ -375,9 +373,6 @@ def mace4D(sino, angles, dist_source_detector, magnification,
end = time.time()
elapsed_t = end-start
print(f"Done computing qGGMRF reconstruction. Elapsed time: {elapsed_t:.2f} sec.")
- if not (save_path is None):
- print("Save qGGMRF reconstruction to disk.")
- np.save(os.path.join(save_path, 'recon_qGGMRF.npy'), init_image)
if np.isscalar(init_image):
[Nz,Nx,Ny], _ = cone3D.compute_img_size(num_views, num_det_rows, num_det_channels,
@@ -401,13 +396,13 @@ def mace4D(sino, angles, dist_source_detector, magnification,
if verbose:
print(f"Begin MACE iteration {itr}/{max_admm_itr}:")
itr_start = time.time()
- # Modify constant_args and variable args respectively for proximal map estimation.
- constant_args['max_iterations'] = max_iterations
- for t in range(Nt):
- variable_args_list[t]['init_image'] = X[0][t]
- variable_args_list[t]['prox_image'] = W[0][t]
# forward model prox map agent
- if not (cluster_ticket is None):
+ if cluster_ticket is not None:
+ # Modify constant_args and variable args respectively for proximal map estimation.
+ constant_args['max_iterations'] = max_iterations
+ for t in range(Nt):
+ variable_args_list[t]['init_image'] = X[0][t]
+ variable_args_list[t]['prox_image'] = W[0][t]
X[0] = np.array(multinode.scatter_gather(cluster_ticket, cone3D.recon,
variable_args_list=variable_args_list,
constant_args=constant_args,
@@ -436,10 +431,6 @@ def mace4D(sino, angles, dist_source_detector, magnification,
denoise_elapsed = denoise_end - denoise_start
print(f"Done denoising in all hyper-planes, elapsed time {denoise_elapsed:.2f} sec")
# save X and W as npy files
- if not (save_path is None):
- for i in range(4):
- np.save(os.path.join(save_path, f'X{i}_itr{itr}.npy'), X[i])
- np.save(os.path.join(save_path, f'W{i}_itr{itr}.npy'), W[i])
Z = sum([beta[k]*(2*X[k]-W[k]) for k in range(4)])
for k in range(4):