tsadakane commented 3 years ago

1. Summary

Backend of CTnoise for Python and addCTnoise for Matlab are replaced with CUDA implementation.

Motivation

I bought new Matlab, but no money to buy Matlab Image Processing Toolbox nor Statictics Toolbox. ToT)

tsadakane commented 3 years ago

2. CUDA implementation

Storategy

Generate poisson and normal random numbers in a kernel.
- Use of curand_poisson and curand_normal of curand.

Verification

Check the distribution of the gerated numbers are the same between CUDA and python.

Code

test_curand.py


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

import _RandomNumberGenerator as RNG

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

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)

nangles = 512 angles = np.linspace(0, 2 * np.pi, nangles, endpoint=False, dtype=np.float32)

fig, axes = plt.subplots(3,2, sharex=True, sharey=True)

list_results = []

####################

GAUSSIAN NOISE

#################### projections_g = np.zeros([nangles, geo.nDetector[0], geo.nDetector[1]], dtype=np.float32) Gaussian = np.array([3, 0.5])

noise_data = np.random.normal(Gaussian[0], Gaussian[1], size=projections_g.shape) # warmup

time_start = time.time() noise_data = np.random.normal(Gaussian[0], Gaussian[1], size=projections_g.shape) time_end = time.time()-time_start print(noise_data.shape) list_results.append("G0 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax = axes[0,0] ax.set_title("Python") ax.set_ylabel("Normal") count, bins, ignored = ax.hist(noise_data, 100, density=True)

time_start = time.time() noise_data = RNG.add_noise(projections_g, Gaussian[0], Gaussian[1]) time_end = time.time()-time_start list_results.append("G1 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax = axes[0,1] ax.set_title("CUDA") count, bins, ignored = ax.hist(noise_data, 100, density=True)

####################

POISSON NOISE

####################

projections_p = np.zeros([nangles, geo.nDetector[0], geo.nDetector[1]], dtype=np.float32) projections_p[:, :, :] = 10 noise_data = np.random.poisson(projections_p) #warmup

time_start = time.time() noise_data = np.random.poisson(projections_p) time_end = time.time()-time_start list_results.append("P0 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax = axes[1, 0] ax.set_ylabel("Poisson") count, bins, ignored = ax.hist(noise_data, 100, density=True)

rng = np.random.default_rng()

time_start = time.time()

noise_data = rng.poisson(projections_p)

time_p[1] = time.time()-time_start

print("P2 mean: {}, sd : {}".format(np.average(noise_data), np.sqrt(np.var(noise_data))))

noise_data = noise_data.reshape([-1])

count, bins, ignored = axes[1, 1].hist(noise_data, 100, density=True)

noise_data = RNG.add_poisson(projections_p) #warmup

time_start = time.time() noise_data = RNG.add_poisson(projections_p) time_end = time.time()-time_start list_results.append("P1 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax=axes[1, 1] count, bins, ignored = ax.hist(noise_data, 100, density=True)

####################

POISSON +GAUSSIAN

####################

noise_data = np.random.poisson(projections_p) # warmup

time_start = time.time() noise_data = np.random.poisson(projections_p) noise_data = noise_data + np.random.normal(Gaussian[0], Gaussian[1], size=noise_data.shape) time_end = time.time()-time_start list_results.append("PG0 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax = axes[2, 0] ax.set_ylabel("Poisson + Normal") count, bins, ignored = ax.hist(noise_data, 100, density=True)

noise_data = RNG.add_noise(projections_p, Gaussian[0], Gaussian[1]) # warmup

time_start = time.time() noise_data = RNG.add_noise(projections_p, Gaussian[0], Gaussian[1]) time_end = time.time()-time_start list_results.append("PG1 time: {}, mean: {}, sd : {}, min: {}, max: {}".format(time_end, np.average(noise_data), np.sqrt(np.var(noise_data)), np.min(noise_data), np.max(noise_data))) noise_data = noise_data.reshape([-1]) ax=axes[2, 1] count, bins, ignored = ax.hist(noise_data, 100, density=True)

####################

Show Results

#################### for line in list_results: print(line)

plt.savefig("CompareCudaAndPython.png") plt.show()

exit()

### Result

* Distributions look the same
* CUDA  version is faster in most of the cases.

G0 time: 0.796917200088501, mean: 3.0000234450711574, sd : 0.4999682128226512, min: 0.3364493374464943, max: 5.64485169303217 G1 time: 0.9529380798339844, mean: 3.0000410079956055, sd : 0.49989402294158936, min: 0.21702957153320312, max: 5.765656471252441 P0 time: 2.7501330375671387, mean: 9.99996104836464, sd : 3.1621095114201148, min: 0, max: 32 P1 time: 0.07698988914489746, mean: 10.00020980834961, sd : 3.1637797355651855, min: 0.0, max: 31.0 PG0 time: 3.5634419918060303, mean: 13.000713100842493, sd : 3.2016939344581368, min: 1.392856437922212, max: 34.41406140912264 PG1 time: 0.0905461311340332, mean: 13.000432968139648, sd : 3.201979160308838, min: 1.1816163063049316, max: 36.43503189086914


![CompareCudaAndPython](https://user-images.githubusercontent.com/40597344/130326552-92e07ad3-87ea-488d-98fe-04687c8e6a7a.png)

### Environment
* Windows 10
* Python 3.8.5 (anaconda)
* MSVC 2019
* RTX 2080Ti x2 + GTX 1070

tsadakane commented 3 years ago

3. Python

Verification

Code

test_NewCTnoise.py

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
from tigre.utilities import CTnoise

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
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

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)

proj_noise = CTnoise.add(proj, Poisson=1e5, Gaussian=np.array([0, 10]))

# Reconstruct
niter = 3
fdk_clean = algs.fdk(proj, geo, angles, gpuids=gpuids)
fdk_noise = algs.fdk(proj_noise, geo, angles, gpuids=gpuids)

# Measure Quality
print("RMSE fdk_clean:")
print(Measure_Quality(fdk_clean, head, ["nRMSE"]))
print("RMSE fdk_noise")
print(Measure_Quality(fdk_noise, head, ["nRMSE"]))

# Plot
fig, axes = plt.subplots(3, 2)
axes[0, 0].set_title("CLEAN")
axes[0, 0].imshow(fdk_clean[geo.nVoxel[0] // 2])
axes[1, 0].imshow(fdk_clean[:, geo.nVoxel[1] // 2, :])
axes[2, 0].imshow(fdk_clean[:, :, geo.nVoxel[2] // 2])
axes[0, 1].set_title("NOISE")
axes[0, 1].imshow(fdk_noise[geo.nVoxel[0] // 2])
axes[1, 1].imshow(fdk_noise[:, geo.nVoxel[1] // 2, :])
axes[2, 1].imshow(fdk_noise[:, :, geo.nVoxel[2] // 2])
plt.savefig("NewCTNoise.png")
plt.show()

Result

Noise is added by CTnoise.add()

RMSE fdk_clean:
429.605506952072
RMSE fdk_noise
525.9876417893452

NewCTNoise_python

Environment

Windows 10
Python 3.8.5 (anaconda)
MSVC 2019
RTX 2080Ti x2 + GTX 1070

tsadakane commented 3 years ago

4. MATLAB

Strategy

As for addCTnoise,...
- keep the default behaviour as it is now.
- use CUDA implementation if the Implementation option is set to cuda
- use CUDA implementation if none of Image Processing Toolbox and Statistics Toolbox are installed even if Implementation is not set.
  Verification
  
  Code 1
  
  Check the distribution

test_curand_mex.m

%% TEST MEX function AddNoise
%% Initialize
clear;
close all;

%% Geometry
listGpuNames = getGpuNames();
gpuids = GpuIds(listGpuNames(1));
geo=defaultGeometry('nVoxel',[256;256;256]);                     
geo.dDetector = [0.8; 0.8] * 2;
geo.sDetector = geo.dDetector.* geo.nDetector;
%disp(geo)

nangles = 512;
angles=linspace(0,2*pi,nangles);

ol_results = [];

% ####################
% ## GAUSSIAN NOISE ##
% ####################
projections_g = zeros(nangles, geo.nDetector(1), geo.nDetector(2), 'single');
Gaussian = [3, 0.5];

% Matlab
time_start = tic;
noise_data = randn(size(projections_g)).* Gaussian(2) + Gaussian(1);
time_end = toc(time_start);
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('G0 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,1);
title('Python');
ylabel('Normal');
histogram(noise_data, 100, 'Normalization', 'probability');

% CUDA
tic;
noise_data = AddNoise(projections_g, Gaussian(1), Gaussian(2));
time_end = toc;
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('G1 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,2);
title('CUDA');
histogram(noise_data, 100, 'Normalization', 'probability');

% ####################
% ## POISSON NOISE  ##
% ####################
projections_p = ones(nangles, geo.nDetector(1), geo.nDetector(2), 'single');
projections_p = projections_p*10;

% Matlab
time_start = tic;
noise_data = poissrnd(projections_p);
time_end = toc(time_start);
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('P0 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,3);
ylabel('Poisson');
histogram(noise_data, 100, 'Normalization', 'probability');

% CUDA
tic;
noise_data = AddNoise(projections_p, 0, 0);
time_end = toc;
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('P1 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,4);
histogram(noise_data, 100, 'Normalization', 'probability');

% ####################
% ## POISSON +GAUSSIAN  ##
% ####################

% Matlab
time_start = tic;
noise_data = poissrnd(projections_p) + randn(size(projections_p)).* Gaussian(2) + Gaussian(1);
time_end = toc(time_start);
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('PG0 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,5);
ylabel('Poisson+Normal');
histogram(noise_data, 100, 'Normalization', 'probability');

% CUDA
tic;
noise_data = AddNoise(projections_p, Gaussian(1), Gaussian(2));
time_end = toc;
noise_data = reshape(noise_data, [1,numel(noise_data)]);
ol_results = [ol_results  sprintf('PG1 time: %f,  mean: %f, sd : %f, min: %f, max: %f\n',time_end, mean(noise_data), sqrt(var(noise_data)), min(noise_data), max(noise_data))];
subplot(3,2,6);
histogram(noise_data, 100, 'Normalization', 'probability');

%
saveas(gcf,'CompareCudaAndMatlab.png')
% Display results
ol_results

Result 1

Distributions look the same
CUDA is much faster especially in generating Poisson distribution random numbers.

G0 time: 0.286028,  mean: 3.000058, sd : 0.500014, min: 0.143546, max: 5.639248
G1 time: 0.125132,  mean: 3.000041, sd : 0.499424, min: 0.217030, max: 5.765656
P0 time: 13.446568,  mean: 9.999731, sd : 3.161444, min: 0.000000, max: 32.000000
P1 time: 0.101630,  mean: 10.000159, sd : 3.163700, min: 0.000000, max: 33.000000
PG0 time: 13.880320,  mean: 12.999625, sd : 3.198814, min: 1.256511, max: 35.258408
PG1 time: 0.089953,  mean: 13.000519, sd : 3.198972, min: 1.181616, max: 36.435032

CompareCudaAndMatlab_R2016a

Environment 1

Windows 10
Matlab R2016a + Image Processing Toolbox + Statistics Toolbox
MSVC 2015
RTX 2080Ti x2 + GTX 1070

Code 2

Check addCTnoise function

test_new_addCTnoise.m

%% TEST New addCTnoise
%% Initialize
clear;
close all;
%% Geometry
listGpuNames = getGpuNames();
gpuids = GpuIds(listGpuNames(1));
geo=defaultGeometry('nVoxel',[128;128;128]);                     

%% Load data and generate projections 
% define angles
nangles = 512;
angles=linspace(0,2*pi,nangles);
% Load thorax phatom data
head=headPhantom(geo.nVoxel);
% generate projections
projections=Ax(head,geo,angles,'interpolated', 'gpuids', gpuids);
% add noise
noise_project_default = addCTnoise(projections);
tic;
noise_projections_cuda  =addCTnoise(projections, 'Implementation', 'cuda');
time_cuda = toc;
tic;
noise_projections_matlab=addCTnoise(projections, 'Implementation', 'matlab');
time_matlab = toc;

%% Reconstruct images using FDK

imgFDK_clean =FDK(projections,geo,angles, 'gpuids', gpuids);
imgFDK_matlab=FDK(noise_projections_matlab,geo,angles, 'gpuids', gpuids);
imgFDK_cuda  =FDK(noise_projections_cuda,geo,angles, 'gpuids', gpuids);

disp(sprintf('TIME  MATLAB: %f sec. CUDA: %f sec.', time_matlab, time_cuda));
disp(sprintf('RMSE  clean: %f  MATLAB: %f, CUDA: %f', ...
    Measure_Quality(imgFDK_clean, head, {'RMSE'}),...
    Measure_Quality(imgFDK_matlab, head, {'RMSE'}),...
    Measure_Quality(imgFDK_cuda, head, {'RMSE'})));

% Show the results
plotImg([imgFDK_clean, imgFDK_matlab,imgFDK_cuda],'Dim','Z', 'Savegif','CompareAddCtNoise.gif');

Result 2

Noise is added by addCTnoise()

TIME  MATLAB: 1.777687 sec. CUDA: 0.277269 sec.
RMSE  clean: 0.024610  MATLAB: 0.026564, CUDA: 0.027060

CompareAddCtNoise_R2016a_Clean_IPToolBox_CUDA

Environment 2

Machine 1
- Windows 10
- Matlab R2016a + Image Processing Toolbox + Statistics Toolbox
- MSVC 2015
- RTX 2080Ti x2 + GTX 1070
machine 2
- Windows 10
- Matlab R2021a + No optional Toolbox
- MSVC2019
- GTX 1060 6GB

tsadakane commented 3 years ago

Discussions

No comparison on the distributions is done between imnoise(poisson) version of Image Processing Toolbox and CUDA version.
- The spec. of the former one is difficult to understand to me.
I seems to me that there is some benefit of changing the default implementation of addCTnoise to cuda, not only from the point of view of time, but also of vanilla-ness.

AnderBiguri commented 3 years ago

Wonderful work, and I completely agree with the vanillaness.

I am away on annual leave the next couple of weeks, I will check this after that and merge!

AnderBiguri commented 3 years ago

Still have not tested this, apologies, will try to do it in the following days. However, I fully agree that making the default CUDA is the best. Feel free to add that.

tsadakane commented 3 years ago

It seems to me that the comments about the default value of Poisson noise (of matlab addCTnoise) needs more explanation.

In order to make this function easier to use, I think that it would be a solution to introduce a parameter that connects the photon count and the ADC of the detector because Poisson noise is based on the photon count and the electric noise is measured in the scale of the ADC. This is another issue though.

AnderBiguri commented 3 years ago

All looks good from here.

Could you do a couple of things?

1- You wrote the code, feel free to remove my name from the coders in the new files 2- Add to Changelog.md the feature description

AnderBiguri commented 3 years ago

In order to make this function easier to use, I think that it would be a solution to introduce a parameter that connects the photon count and the ADC of the detector because Poisson noise is based on the photon count and the electric noise is measured in the scale of the ADC. This is another issue though.

I do agree here, but I never knew how to add this. In my experience, the detectors (and users) tend to scale the max photon count to the max value of the ADC, so there is no direct way to know the photon count from a measured projection set. I agree that the description of I0 should be better though.

tsadakane commented 3 years ago

so there is no direct way to know the photon count from a measured projection set

I agree. So I think the parameter value will come from outside. The detector manufacturer knows [ADC Counts/uR] for the specific mode of the specific detector. Many radiologists measure the value (and the spectrum) by themselves (I guess). If the projection data is based on computer simulation, the user can choose the value (the user might be able to add noise by him(her)self instead of using TIGRE in this case though.). As for the photon counting detector, I think that the projection data is nothing but the photon count, but I am not sure about the electric noise in this case...

So far, I don't have an idea that is compatible with existing code though... :-P

AnderBiguri commented 3 years ago

I suggest we leave I0 as is then, as I don't see a way to cater to all different detector types...

When I Worked with uCT scans, there in each single scan, the photon count was different, for example. You could analogically tune kV and mA for each scan, so there wasn't even a possibility of knowing how much there was....

tsadakane commented 3 years ago

I suggest we leave I0 as is then, as I don't see a way to cater to all different detector types...

Agree. I am not going add commit for this PR.

You could analogically tune kV and mA for each scan, so there wasn't even a possibility of knowing how much there was....

True. So it is necessary to do experiments using a device like https://www.amptek.com/products/x-ray-detectors/cdte-x-ray-and-gamma-ray-detectors/cdte-x-ray-and-gamma-ray-detector. The output of this device is the numbers of the photons of each energy bins, detected in the certain time window.

tsadakane commented 3 years ago

Thank you for accepting my PR and the discussion.

AnderBiguri commented 3 years ago

Absolutely! I really appreciate all the fantastic work you are contributing, I can't thank you enough :)

CERN / TIGRE

CUDA implementation of CTnoise #320

1. Summary

Motivation

2. CUDA implementation

Storategy

Verification

Code

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.

Geometry

print(geo)

GAUSSIAN NOISE

POISSON NOISE

rng = np.random.default_rng()

time_start = time.time()

noise_data = rng.poisson(projections_p)

time_p[1] = time.time()-time_start

print("P2 mean: {}, sd : {}".format(np.average(noise_data), np.sqrt(np.var(noise_data))))

noise_data = noise_data.reshape([-1])

count, bins, ignored = axes[1, 1].hist(noise_data, 100, density=True)

POISSON +GAUSSIAN

Show Results

3. Python

Verification

Code

Result

Environment

4. MATLAB

Strategy

Verification

Code 1

Result 1

Environment 1

Code 2

Result 2

Environment 2