From bbb98660003f65a31dd464f11286ddf842311b2f Mon Sep 17 00:00:00 2001
From: Charles Bouman <Charles.Bouman@gmail.com>
Date: Sat, 4 Jun 2022 18:34:03 -0400
Subject: [PATCH 1/4] Changed default iterations to 100
---
demo/demo_3D_shepp_logan.py | 41 ++++++++++++++++++++++---------------
mbircone/cone3D.py | 4 ++--
2 files changed, 27 insertions(+), 18 deletions(-)
diff --git a/demo/demo_3D_shepp_logan.py b/demo/demo_3D_shepp_logan.py
index 61779ec..c5be708 100644
--- a/demo/demo_3D_shepp_logan.py
+++ b/demo/demo_3D_shepp_logan.py
@@ -31,9 +31,15 @@ num_views = 64
# projection angles will be uniformly spaced within the range [0, 2*pi).
angles = np.linspace(0, 2 * np.pi, num_views, endpoint=False)
# qGGMRF recon parameters
-sharpness = 2.0 # Controls regularization level of reconstruction
-max_iterations=40
-stop_threshold=0.005
+sharpness = 1.0 # Controls regularization level of reconstruction by controling prior term weighting
+snr_db = 40 # Controls regularization level of reconstruction by controling data term weighting
+T = 0.1 # Controls edginess of reconstruction
+# convergence parameters
+stop_threshold = 0.005
+# display parameters
+vmin = 0.10
+vmax = 0.12
+
# local path to save phantom, sinogram, and reconstruction images
save_path = f'output/3D_shepp_logan/'
os.makedirs(save_path, exist_ok=True)
@@ -57,8 +63,10 @@ print('ROI shape is:', num_slices_ROI, num_rows_ROI, num_cols_ROI)
######################################################################################
# Generate a 3D shepp logan phantom with size calculated in previous section
######################################################################################
-phantom = mbircone.phantom.gen_shepp_logan_3d(num_rows_ROI, num_cols_ROI, num_slices_ROI) # generate phantom within ROI
+phantom = mbircone.phantom.gen_shepp_logan_3d(num_rows_ROI, num_cols_ROI, num_slices_ROI) # generate phantom within ROI
phantom = mbircone.cone3D.pad_roi2ror(phantom, boundary_size) # zero-pad phantom to ROR
+# scale the phantom by a factor of 10.0 to make the projections physical realistic -log attenuation values
+phantom = phantom/10.0
print('Padded ROR phantom shape = ', np.shape(phantom))
@@ -76,17 +84,14 @@ print('Synthetic sinogram shape: (num_views, num_det_rows, num_det_channels) = '
# Perform 3D qGGMRF reconstruction
######################################################################################
print('Performing 3D qGGMRF reconstruction ...')
-recon = mbircone.cone3D.recon(sino, angles,
- dist_source_detector, magnification,
- sharpness=sharpness,
- max_iterations=max_iterations, stop_threshold=stop_threshold)
+recon = mbircone.cone3D.recon(sino, angles, dist_source_detector, magnification, sharpness=sharpness, snr_db = snr_db, T = T, stop_threshold = stop_threshold)
print('recon shape = ', np.shape(recon))
######################################################################################
# Generate phantom, synthetic sinogram, and reconstruction images
######################################################################################
-# Generate sino and recon images
+# sinogram images
for view_idx in [0, num_views//4, num_views//2]:
view_angle = int(angles[view_idx]*180/np.pi)
plot_image(sino[view_idx, :, :], title=f'sinogram view angle {view_angle} ',
@@ -97,17 +102,21 @@ display_x = num_rows_ROR // 2
display_y = num_cols_ROR // 2
# phantom images
plot_image(phantom[display_slice], title=f'phantom, axial slice {display_slice}',
- filename=os.path.join(save_path, 'phantom_axial.png'), vmin=1.0, vmax=1.1)
-plot_image(phantom[:,display_x,:], title=f'phantom, coronal slice {display_x}',
- filename=os.path.join(save_path, 'phantom_coronal.png'), vmin=1.0, vmax=1.1)
+ filename=os.path.join(save_path, 'phantom_axial.png'), vmin=0, vmax=0.40)
+plot_image(phantom[display_slice], title=f'phantom, axial slice {display_slice}',
+ filename=os.path.join(save_path, 'phantom_axial.png'), vmin=vmin, vmax=vmax)
+plot_image(phantom[:,display_x,:], title=f'phantom, coronal slice {display_x}',
+ filename=os.path.join(save_path, 'phantom_coronal.png'), vmin=vmin, vmax=vmax)
plot_image(phantom[:,:,display_y], title=f'phantom, sagittal slice {display_y}',
- filename=os.path.join(save_path, 'phantom_sagittal.png'), vmin=1.0, vmax=1.1)
+ filename=os.path.join(save_path, 'phantom_sagittal.png'), vmin=vmin, vmax=vmax)
# recon images
plot_image(recon[display_slice], title=f'qGGMRF recon, axial slice {display_slice}',
- filename=os.path.join(save_path, 'recon_axial.png'), vmin=1.0, vmax=1.1)
+ filename=os.path.join(save_path, 'recon_axial.png'), vmin=0, vmax=0.40)
+plot_image(recon[display_slice], title=f'qGGMRF recon, axial slice {display_slice}',
+ filename=os.path.join(save_path, 'recon_axial.png'), vmin=vmin, vmax=vmax)
plot_image(recon[:,display_x,:], title=f'qGGMRF recon, coronal slice {display_x}',
- filename=os.path.join(save_path, 'recon_coronal.png'), vmin=1.0, vmax=1.1)
+ filename=os.path.join(save_path, 'recon_coronal.png'), vmin=vmin, vmax=vmax)
plot_image(recon[:,:,display_y], title=f'qGGMRF recon, sagittal slice {display_y}',
- filename=os.path.join(save_path, 'recon_sagittal.png'), vmin=1.0, vmax=1.1)
+ filename=os.path.join(save_path, 'recon_sagittal.png'), vmin=vmin, vmax=vmax)
print(f"Images saved to {save_path}.")
input("Press Enter")
diff --git a/mbircone/cone3D.py b/mbircone/cone3D.py
index 4d182f7..e05e1b1 100644
--- a/mbircone/cone3D.py
+++ b/mbircone/cone3D.py
@@ -486,7 +486,7 @@ def recon(sino, angles, dist_source_detector, magnification,
init_image=0.0, prox_image=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=20, stop_threshold=0.02,
+ sharpness=0.0, sigma_x=None, sigma_p=None, max_iterations=100, stop_threshold=0.02,
num_threads=None, NHICD=False, verbose=1, lib_path=__lib_path):
"""Compute 3D cone beam MBIR reconstruction
@@ -541,7 +541,7 @@ def recon(sino, angles, dist_source_detector, magnification,
sigma_p (float, optional): [Default=None] Scalar value :math:`>0` that specifies the proximal map parameter.
Ignored if prox_image is None.
If None and proximal image is not None, automatically set with auto_sigma_p. Regularization should be controled with the ``sharpness`` parameter, but ``sigma_p`` can be set directly by expert users.
- max_iterations (int, optional): [Default=20] Integer valued specifying the maximum number of iterations.
+ max_iterations (int, optional): [Default=100] Integer valued specifying the maximum number of iterations.
stop_threshold (float, optional): [Default=0.02] Scalar valued stopping threshold in percent.
If stop_threshold=0.0, then run max iterations.
num_threads (int, optional): [Default=None] Number of compute threads requested when executed.
--
GitLab
From ef15bb41490a0005c0a45fec36545f9c4d8c1f3e Mon Sep 17 00:00:00 2001
From: dyang37 <yang1467@purdue.edu>
Date: Sun, 5 Jun 2022 14:48:55 -0400
Subject: [PATCH 2/4] change default snr_db value from 30.0 to 40.0
---
mbircone/cone3D.py | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/mbircone/cone3D.py b/mbircone/cone3D.py
index e05e1b1..1e09e12 100644
--- a/mbircone/cone3D.py
+++ b/mbircone/cone3D.py
@@ -122,7 +122,7 @@ def calc_weights(sino, weight_type):
return weights
-def auto_sigma_y(sino, magnification, weights, snr_db=30.0, delta_pixel_image=1.0, delta_pixel_detector=1.0):
+def auto_sigma_y(sino, magnification, weights, snr_db=40.0, delta_pixel_image=1.0, delta_pixel_detector=1.0):
"""Compute the automatic value of ``sigma_y`` for use in MBIR reconstruction.
Args:
@@ -132,7 +132,7 @@ def auto_sigma_y(sino, magnification, weights, snr_db=30.0, delta_pixel_image=1.
numpy array of weights with same shape as sino.
The parameters weights should be the same values as used in mbircone reconstruction.
snr_db (float, optional):
- [Default=30.0] Scalar value that controls assumed signal-to-noise ratio of the data in dB.
+ [Default=40.0] Scalar value that controls assumed signal-to-noise ratio of the data in dB.
delta_pixel_image (float, optional):
[Default=1.0] Scalar value of pixel spacing in :math:`ALU`.
delta_pixel_detector (float, optional):
@@ -484,7 +484,7 @@ def recon(sino, angles, dist_source_detector, magnification,
channel_offset=0.0, row_offset=0.0, rotation_offset=0.0,
delta_pixel_detector=1.0, delta_pixel_image=None, ror_radius=None,
init_image=0.0, prox_image=None,
- sigma_y=None, snr_db=30.0, weights=None, weight_type='unweighted',
+ sigma_y=None, snr_db=40.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=100, stop_threshold=0.02,
num_threads=None, NHICD=False, verbose=1, lib_path=__lib_path):
@@ -513,7 +513,7 @@ def recon(sino, angles, dist_source_detector, magnification,
sigma_y (float, optional): [Default=None] Scalar value of noise standard deviation parameter.
If None, automatically set with auto_sigma_y.
- snr_db (float, optional): [Default=30.0] Scalar value that controls assumed signal-to-noise ratio of the data in dB.
+ snr_db (float, optional): [Default=40.0] Scalar value that controls assumed signal-to-noise ratio of the data in dB.
Ignored if sigma_y is not None.
weights (ndarray, optional): [Default=None] 3D weights array with same shape as sino.
weight_type (string, optional): [Default='unweighted'] Type of noise model used for data.
--
GitLab
From 3c8303c966cd1f087630dc9a0c94dbcada4df69d Mon Sep 17 00:00:00 2001
From: dyang37 <yang1467@purdue.edu>
Date: Sun, 5 Jun 2022 14:50:20 -0400
Subject: [PATCH 3/4] remove snr_db variable from 3D shepp logan demo, since
now the default value is 40.0.
---
demo/demo_3D_shepp_logan.py | 3 +--
1 file changed, 1 insertion(+), 2 deletions(-)
diff --git a/demo/demo_3D_shepp_logan.py b/demo/demo_3D_shepp_logan.py
index c5be708..d361d4c 100644
--- a/demo/demo_3D_shepp_logan.py
+++ b/demo/demo_3D_shepp_logan.py
@@ -32,7 +32,6 @@ num_views = 64
angles = np.linspace(0, 2 * np.pi, num_views, endpoint=False)
# qGGMRF recon parameters
sharpness = 1.0 # Controls regularization level of reconstruction by controling prior term weighting
-snr_db = 40 # Controls regularization level of reconstruction by controling data term weighting
T = 0.1 # Controls edginess of reconstruction
# convergence parameters
stop_threshold = 0.005
@@ -84,7 +83,7 @@ print('Synthetic sinogram shape: (num_views, num_det_rows, num_det_channels) = '
# Perform 3D qGGMRF reconstruction
######################################################################################
print('Performing 3D qGGMRF reconstruction ...')
-recon = mbircone.cone3D.recon(sino, angles, dist_source_detector, magnification, sharpness=sharpness, snr_db = snr_db, T = T, stop_threshold = stop_threshold)
+recon = mbircone.cone3D.recon(sino, angles, dist_source_detector, magnification, sharpness=sharpness, T = T, stop_threshold = stop_threshold)
print('recon shape = ', np.shape(recon))
--
GitLab
From 57f181641b49d87ec9204d81e15057438293e8dc Mon Sep 17 00:00:00 2001
From: Charles Bouman <Charles.Bouman@gmail.com>
Date: Sun, 5 Jun 2022 16:29:40 -0400
Subject: [PATCH 4/4] Reduce default blocksize parameter for shepp-logan
phantom
---
mbircone/phantom.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/mbircone/phantom.py b/mbircone/phantom.py
index 4213247..757b97d 100644
--- a/mbircone/phantom.py
+++ b/mbircone/phantom.py
@@ -77,7 +77,7 @@ def gen_microscopy_sample(num_rows, num_cols):
return image
-def gen_shepp_logan_3d(num_rows, num_cols, num_slices, block_size=(4,4,4)):
+def gen_shepp_logan_3d(num_rows, num_cols, num_slices, block_size=(2,2,2)):
"""
Generate a smoothed 3D Shepp Logan phantom with block-averaging strategy.
--
GitLab