Enh: Python order subset

[tsadakane]

A solution for #268. #313 is related.

• [x] Add plot_angle function
• [x] Modify order_subsets to return lists of np.ndarray's
• [x] Modify os-xxxx

[tsadakane]

Add plot_angle function.

The following program works with the parameters written in it. It may not work as expected with different parameters, strategy="ramdon" for example. This is because order_subsets returns different type values depending on the mode paramter. This will be fixed later.

#%%
import numpy as np
import importlib
import tigre.utilities.order_subsets
from tigre.utilities.visualization.plot_angles import plot_angles

#%%
importlib.reload(tigre.utilities.order_subsets)
# %%
nangles = 16
range_over_pi = 2
angles = np.linspace(0, range_over_pi * np.pi, nangles, endpoint=False, dtype=np.float32)
# if angles.ndim == 1:
#     a1 = angles
#     a2 = np.zeros(angles.shape[0], dtype=np.float32)
#     angles = np.vstack((a1, a2, a2)).T
print(f"angles.shape = {angles.shape}")
print(f"angles.ndim  = {angles.ndim}")

# %%
block_size = 1
strategy = "angularDistance"# "random2"#"random" #"ordered" #"angularDistance"
print(f"block_size = {block_size}")
print(f"strategy   = {strategy}")
print(f"nangles    = {nangles}")
print(f"range      = {range_over_pi} PI")

angle_blocks, angle_index = tigre.utilities.order_subsets.order_subsets(angles, blocksize= block_size, mode=strategy)

# %%
filename = f"ordered_subsets_{range_over_pi}PI_nangles_{nangles}_blocksize_{block_size}_mode_{strategy}.gif"
plot_angles(angle_blocks, angle_index, savegif=filename, angles_extended = (angles.ndim != 1))

[tsadakane]

Modify order_subsets

#%%
import numpy as np
import importlib
import tigre.utilities.order_subsets
from tigre.utilities.visualization.plot_angles import plot_angles

#%%
importlib.reload(tigre.utilities.order_subsets)
# %%
nangles = 17
range_over_pi = 2
block_size = 4
strategy = "ordered"# "random2"#"random" #"ordered" #"angularDistance"

#%%
angles = np.linspace(0, range_over_pi * np.pi, nangles, endpoint=False, dtype=np.float32)
if angles.ndim == 1:
    a1 = angles
    a2 = np.zeros(angles.shape[0], dtype=np.float32)
    angles = np.vstack((a1, a2, a2)).T
print(f"angles.shape = {angles.shape}")
print(f"angles.ndim  = {angles.ndim}")
print(f"block_size   = {block_size}")
print(f"strategy     = {strategy}")
print(f"nangles      = {nangles}")
print(f"range        = {range_over_pi} PI")

angle_blocks, angle_index = tigre.utilities.order_subsets.order_subsets(angles, blocksize= block_size, mode=strategy)

# %%
filename = f"ordered_subsets_{range_over_pi}PI_nangles_{nangles}_blocksize_{block_size}_mode_{strategy}.gif"
plot_angles(angle_blocks, angle_index, savegif=filename, angles_extended = (angles.ndim != 1))

The code above generates the following figure:

Modifying as

nangles = 17
range_over_pi = 1
block_size = 2
strategy = "angularDistance"# "random2"#"random" #"ordered" #"angularDistance"

it generates

[tsadakane]

Modify os-xxx

The following code is a modification of example.py. Instead of FDK, OSSART with blocksize = 17 is used.

import numpy as np
import tigre
import tigre.algorithms as algs
from tigre.utilities import sample_loader
from tigre.utilities.Measure_Quality import Measure_Quality
import tigre.utilities.gpu as gpu
import matplotlib.pyplot as plt

### This is just a basic example of very few TIGRE functionallity.
# We hihgly recomend checking the Demos folder, where most if not all features of tigre are demoed.

listGpuNames = gpu.getGpuNames()
if len(listGpuNames) == 0:
    print("Error: No gpu found")
else:
    for id in range(len(listGpuNames)):
        print("{}: {}".format(id, listGpuNames[id]))

gpuids = gpu.getGpuIds(listGpuNames[0])
print(gpuids)

# Geometry
# geo1 = tigre.geometry(mode='cone', high_resolution=False, default=True)
geo = tigre.geometry(mode="cone", nVoxel=np.array([256, 256, 256]), default=True)
geo.dDetector = np.array([0.8, 0.8]) * 2  # size of each pixel            (mm)
geo.sDetector = geo.dDetector * geo.nDetector
# print(geo)

# blocksize = 17
nangles = 100
angles = np.linspace(0, 2 * np.pi, nangles, endpoint=False, dtype=np.float32)

# Prepare projection data
head = sample_loader.load_head_phantom(geo.nVoxel)
proj = tigre.Ax(head, geo, angles, gpuids=gpuids)

# Reconstruct
niter = 20 #5
out1 = algs.ossart(proj, geo, angles, niter, blocksize=17, gpuids=gpuids)
out2 = algs.ossart(proj, geo, angles, niter, blocksize=20, gpuids=gpuids)

# Measure Quality
# 'RMSE', 'MSSIM', 'SSD', 'UQI'
print("RMSE ossart 1:")
print(Measure_Quality(out1, head, ["nRMSE"]))
print("RMSE ossart 2")
print(Measure_Quality(out2, head, ["nRMSE"]))

# Plot
fig, axes = plt.subplots(3, 2)
axes[0, 0].set_title("OS-SART 1")
axes[0, 0].imshow(out1[geo.nVoxel[0] // 2])
axes[1, 0].imshow(out1[:, geo.nVoxel[1] // 2, :])
axes[2, 0].imshow(out1[:, :, geo.nVoxel[2] // 2])
axes[0, 1].set_title("OS-SART 2")
axes[0, 1].imshow(out2[geo.nVoxel[0] // 2])
axes[1, 1].imshow(out2[:, geo.nVoxel[1] // 2, :])
axes[2, 1].imshow(out2[:, :, geo.nVoxel[2] // 2])
plt.show()
# tigre.plotProj(proj)
# tigre.plotImg(out1)

It displayed on the console

>python testOSSART.py
0: NVIDIA GeForce GTX 1060 6GB
{'name': 'NVIDIA GeForce GTX 1060 6GB', 'devices': [0]}
order_subsets(blocksize = 17, mode = None)  <<<<< print is added at the head of order_subsets
OSSART algorithm in progress.
Estimated time until completion : 00:00:32
order_subsets(blocksize = 20, mode = None)  <<<<< print is added at the head of order_subsets
OSSART algorithm in progress.
Estimated time until completion : 00:00:28
RMSE ossart 1:
236.7852370511056
RMSE ossart 2
264.2052417134304

and the reconstructed images were

[AnderBiguri]

Great!

I assume random2 is the random ordering, but where the block angles themselves are also random, instead of contiguous, right? Good to have it!

Can you, in this PR, also move all the visualization functions to the new folder?

[tsadakane]

Thank you for your comment.

Can you, in this PR, also move all the visualization functions to the new folder?

Yes, I can move plotImg, plotProj, plot_geometry into utilities/visualizaion, but most of the demo scripts must be modified at the same time and it is not related to order_subset issue. I think it would be better to do it in another PR. What do you think?

[AnderBiguri]

Ah yes! you are correct. Sorry I had my MATLAB brain on, where no modification would be needed. But indeed, in python, they would. So I agree.