From 414ca427b0206396aa725524a66d295b81b844cc Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Sat, 21 Aug 2021 22:58:13 +0900
Subject: [PATCH 1/7] Add CUDA implementation of random number generator

* Added Common files
---
 Common/CUDA/RandomNumberGenerator.cu  | 194 ++++++++++++++++++++++++++
 Common/CUDA/RandomNumberGenerator.hpp |  50 +++++++
 2 files changed, 244 insertions(+)
 create mode 100644 Common/CUDA/RandomNumberGenerator.cu
 create mode 100644 Common/CUDA/RandomNumberGenerator.hpp

diff --git a/Common/CUDA/RandomNumberGenerator.cu b/Common/CUDA/RandomNumberGenerator.cu
new file mode 100644
index 00000000..f9cd2188
--- /dev/null
+++ b/Common/CUDA/RandomNumberGenerator.cu
@@ -0,0 +1,194 @@
+/*-------------------------------------------------------------------------
+ *
+ * CUDA functions for random number generator
+ *
+ * Adds noise of Poisson and normal distribution to the input.
+ *
+ * CODE by       Ander Biguri
+ *               Tomoyuki SADAKANE
+ * ---------------------------------------------------------------------------
+ * ---------------------------------------------------------------------------
+ * Copyright (c) 2015, University of Bath and CERN- European Organization for
+ * Nuclear Research
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ *
+ * 3. Neither the name of the copyright holder nor the names of its contributors
+ * may be used to endorse or promote products derived from this software without
+ * specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ * ---------------------------------------------------------------------------
+ *
+ * Contact: tigre.toolbox@gmail.com
+ * Codes  : https://github.com/CERN/TIGRE
+ * ---------------------------------------------------------------------------
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <cuda.h>
+#include <curand_kernel.h>
+#include <curand.h>
+
+#include "gpuUtils.hpp"
+#include "RandomNumberGenerator.hpp"
+
+#define cudaCheckErrors(msg) \
+do { \
+        cudaError_t __err = cudaGetLastError(); \
+        if (__err != cudaSuccess) { \
+                mexPrintf("%s \n",msg);\
+                cudaDeviceReset();\
+                mexErrMsgIdAndTxt("RandomNumberGenerator:",cudaGetErrorString(__err));\
+        } \
+} while (0)
+
+
+__global__ void setup_kernel(curandState *state) {
+    int idx = threadIdx.x + blockIdx.x * blockDim.x;
+    /* Each thread gets same seed, a different sequence number, no offset */
+    curand_init(1234, idx, 0, &state[idx]);
+}
+
+__global__ void GeneratePoisson(curandState *state, const float* pfIn, size_t uiLen, float* pfOut) {
+    int idx = threadIdx.x + blockIdx.x * blockDim.x;
+    /* Copy state to local memory for efficiency */
+    curandState localState = state[idx];
+    int iIter = (uiLen + blockDim.x*gridDim.x - 1)/(blockDim.x*gridDim.x);
+    for (int iI = 0; iI < iIter; ++iI) {
+        size_t uiPos = (size_t)blockDim.x*gridDim.x*iI+idx;
+        if (uiPos < uiLen) {
+            /* Poisson */
+            unsigned int uiPoisson = curand_poisson(&localState, pfIn[uiPos]);
+            pfOut[uiPos] = (float)uiPoisson;
+        }
+    }
+    /* Copy state back to global memory */
+    state[idx] = localState;
+}
+
+__global__ void GeneratePoissonAddGaussian(curandState *state,
+                        const float* pfIn,
+                        size_t uiLen, 
+                        float fGaussMu,
+                        float fGaussSigma,
+                        float* pfOut)
+{
+    int idx = threadIdx.x + blockIdx.x * blockDim.x;
+    /* Copy state to local memory for efficiency */
+    curandState localState = state[idx];
+    int iIter = (uiLen + blockDim.x*gridDim.x - 1)/(blockDim.x*gridDim.x);
+    for (int iI = 0; iI < iIter; ++iI) {
+        size_t uiPos = (size_t)blockDim.x*gridDim.x*iI+idx;
+        if (uiPos < uiLen) {
+            /* Poisson */
+            unsigned int uiPoisson = curand_poisson(&localState, pfIn[uiPos]);
+            /* Gaussian */
+            float fNormal = curand_normal(&localState) * fGaussSigma + fGaussMu;
+            pfOut[uiPos] = fNormal + (float)uiPoisson;
+        }
+    }
+    /* Copy state back to global memory */
+    state[idx] = localState;
+}
+
+
+template<class T_value>
+void GetMinMax(const T_value* pfIn, size_t uiLen, T_value& tvMin, T_value& tvMax) {
+    tvMin = pfIn[0];
+    tvMax = pfIn[0];
+    T_value tvVal;
+    for (int iI = 1; iI < uiLen; ++iI) {
+        tvVal = pfIn[iI];
+        if (tvMax < tvVal) { tvMax = tvVal; continue;}
+        if (tvMin > tvVal) { tvMin = tvVal; continue;}
+    }
+}
+void poisson_1d(const float* pfIn, size_t uiLen, float* pfOut, const GpuIds& gpuids) {
+    // printf("poisson_1d(pfIn = %p, uiLen = %zd, pfOut = %p)\n", pfIn, uiLen, pfOut);
+    float* d_pfIn = nullptr;
+    float* d_pfOut = nullptr;
+    cudaMalloc((void **)&d_pfIn, uiLen * sizeof(float));
+    cudaCheckErrors("poisson_1d fail cudaMalloc 1");
+    cudaMalloc((void **)&d_pfOut, uiLen * sizeof(float));
+    cudaCheckErrors("poisson_1d fail cudaMalloc 2");
+    cudaMemcpy(d_pfIn, pfIn, uiLen*sizeof(float), cudaMemcpyHostToDevice);
+    cudaCheckErrors("poisson_1d fail cudaMemcpy 1");
+
+    // float fMin, fMax;
+    // GetMinMax(pfIn, uiLen, fMin, fMax);
+    // printf("fMin, fMax = %f, %f\n", fMin, fMax);
+    curandState *curandStates = nullptr;
+    const int kiBlockDim = 1024;  // Threads per Block
+    const int kiGridDim = 64;//(uiLen+kiBlockDim-1)/kiBlockDim;
+    cudaMalloc((void **)&curandStates, kiGridDim * kiBlockDim * sizeof(curandState));
+    cudaCheckErrors("poisson_1d fail cudaMalloc 3");
+    setup_kernel<<<kiGridDim, kiBlockDim>>>(curandStates);
+    GeneratePoisson<<<kiGridDim, kiBlockDim>>>(curandStates, d_pfIn, uiLen, d_pfOut);
+    cudaMemcpy(pfOut, d_pfOut, uiLen*sizeof(float), cudaMemcpyDeviceToHost);
+    cudaCheckErrors("poisson_1d fail cudaMemcpy 2");
+    // GetMinMax(pfOut, uiLen, fMin, fMax);
+    // printf("fMin, fMax = %f, %f\n", fMin, fMax);
+    
+    cudaFree(d_pfIn); d_pfIn = nullptr;
+    cudaFree(d_pfOut); d_pfOut = nullptr;
+    cudaFree(curandStates); curandStates = nullptr;
+}
+
+void poisson_gaussian_1d(const float* pfIn,
+                        size_t uiLen,
+                        float fGaussMu,
+                        float fGaussSigma,
+                        float* pfOut,
+                        GpuIds& gpuids)
+{
+    // printf("poisson_gaussian_1d(pfIn = %p, uiLen = %zd, fGaussMu = %+f, fGaussSigma = %f, pfOut = %p)\n", pfIn, uiLen, fGaussMu, fGaussSigma, pfOut);
+    float* d_pfIn = nullptr;
+    float* d_pfOut = nullptr;
+    cudaMalloc((void **)&d_pfIn, uiLen * sizeof(float));
+    cudaCheckErrors("poisson_gaussian_1d fail cudaMalloc 1");
+    cudaMalloc((void **)&d_pfOut, uiLen * sizeof(float));
+    cudaCheckErrors("poisson_gaussian_1d fail cudaMalloc 2");
+    cudaMemcpy(d_pfIn, pfIn, uiLen*sizeof(float), cudaMemcpyHostToDevice);
+    cudaCheckErrors("poisson_gaussian_1d fail cudaMemcpy 1");
+
+    // float fMin, fMax;
+    // GetMinMax(pfIn, uiLen, fMin, fMax);
+    // printf("fMin, fMax = %f, %f\n", fMin, fMax);
+    curandState *curandStates = nullptr;
+    const int kiBlockDim = 64;  // Threads per Block
+    const int kiGridDim = 64;//(uiLen+kiBlockDim-1)/kiBlockDim;
+    cudaMalloc((void **)&curandStates, kiGridDim * kiBlockDim * sizeof(curandState));
+    cudaCheckErrors("poisson_gaussian_1d fail cudaMalloc 3");
+    setup_kernel<<<kiGridDim, kiBlockDim>>>(curandStates);
+    GeneratePoissonAddGaussian<<<kiGridDim, kiBlockDim>>>(curandStates, d_pfIn, uiLen, fGaussMu, fGaussSigma, d_pfOut);
+    cudaMemcpy(pfOut, d_pfOut, uiLen*sizeof(float), cudaMemcpyDeviceToHost);
+    cudaCheckErrors("poisson_gaussian_1d fail cudaMemcpy 2");
+    // GetMinMax(pfOut, uiLen, fMin, fMax);
+    // printf("fMin, fMax = %f, %f\n", fMin, fMax);
+
+
+    cudaFree(d_pfIn); d_pfIn = nullptr;
+    cudaFree(d_pfOut); d_pfOut = nullptr;
+    cudaFree(curandStates); curandStates = nullptr;
+}
diff --git a/Common/CUDA/RandomNumberGenerator.hpp b/Common/CUDA/RandomNumberGenerator.hpp
new file mode 100644
index 00000000..d6f6c6c2
--- /dev/null
+++ b/Common/CUDA/RandomNumberGenerator.hpp
@@ -0,0 +1,50 @@
+/*-------------------------------------------------------------------------
+ *
+ * Header CUDA functions for random number generator
+ *
+ * Adds noise of Poisson and normal distribution to the input.
+ *
+ * CODE by       Ander Biguri
+ *               Tomoyuki SADAKANE
+ * ---------------------------------------------------------------------------
+ * ---------------------------------------------------------------------------
+ * Copyright (c) 2015, University of Bath and CERN- European Organization for
+ * Nuclear Research
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice,
+ * this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ * this list of conditions and the following disclaimer in the documentation
+ * and/or other materials provided with the distribution.
+ *
+ * 3. Neither the name of the copyright holder nor the names of its contributors
+ * may be used to endorse or promote products derived from this software without
+ * specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ * ---------------------------------------------------------------------------
+ *
+ * Contact: tigre.toolbox@gmail.com
+ * Codes  : https://github.com/CERN/TIGRE
+ * ---------------------------------------------------------------------------
+ */
+
+#include "TIGRE_common.hpp"
+#include "GpuIds.hpp"
+void poisson_1d(const float* pfIn, size_t uiLen, float* pfOut, const GpuIds& gpuids);
+void poisson_gaussian_1d(const float* pfPoissonL, size_t uiLen, float fGaussMu, float fGaussSigma, float* pfOut, GpuIds& gpuids);

From 1b4bb498cfdfc34a581724fbc4ee182aedce5806 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Sat, 21 Aug 2021 23:38:49 +0900
Subject: [PATCH 2/7] Add CUDA implementation of random number generator

* Modified and added Matlab files
---
 MATLAB/Compile.m                             |   6 +
 MATLAB/Utilities/addCTnoise.m                |  46 ++++---
 MATLAB/Utilities/cuda_interface/AddNoise.cpp | 126 +++++++++++++++++++
 3 files changed, 161 insertions(+), 17 deletions(-)
 create mode 100644 MATLAB/Utilities/cuda_interface/AddNoise.cpp

diff --git a/MATLAB/Compile.m b/MATLAB/Compile.m
index 0390b315..c2d3054e 100644
--- a/MATLAB/Compile.m
+++ b/MATLAB/Compile.m
@@ -61,6 +61,7 @@
         mex -largeArrayDims ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
         mex -largeArrayDims ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
         mex -largeArrayDims ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
+        mex -largeArrayDims ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
         mex -largeArrayDims ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/win64
         mex -largeArrayDims ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
         mex -largeArrayDims ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win64
@@ -70,6 +71,7 @@
         mex  ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
         mex  ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
         mex  ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
+        mex  ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
         mex  ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/win32
         mex  ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
         mex  ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/win32
@@ -83,6 +85,7 @@
         mex -largeArrayDims ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
         mex -largeArrayDims ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
         mex -largeArrayDims ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
+        mex -largeArrayDims ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
         mex -largeArrayDims ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/mac64
         mex -largeArrayDims ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
         mex -largeArrayDims ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac64
@@ -92,6 +95,7 @@
         mex  ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
         mex  ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
         mex  ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
+        mex  ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
         mex  ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/mac32
         mex  ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
         mex  ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/mac32
@@ -104,6 +108,7 @@
         mex -largeArrayDims ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
         mex -largeArrayDims ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
         mex -largeArrayDims ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
+        mex -largeArrayDims ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
         mex -largeArrayDims ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/linux64
         mex -largeArrayDims ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
         mex -largeArrayDims ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux64
@@ -113,6 +118,7 @@
         mex  ./Utilities/cuda_interface/minTV.cpp ../Common/CUDA/POCS_TV.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
         mex  ./Utilities/cuda_interface/AwminTV.cpp ../Common/CUDA/POCS_TV2.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
         mex  ./Utilities/cuda_interface/tvDenoise.cpp ../Common/CUDA/tvdenoising.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
+        mex  ./Utilities/cuda_interface/AddNoise.cpp ../Common/CUDA/RandomNumberGenerator.cu ../Common/CUDA/GpuIds.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
         mex  -largeArrayDims ./Utilities/IO/VarianCBCT/mexReadXim.cpp -outdir ./Mex_files/linux32
         mex  ./Utilities/GPU/getGpuName_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
         mex  ./Utilities/GPU/getGpuCount_mex.cpp ../Common/CUDA/gpuUtils.cu -outdir ./Mex_files/linux32
diff --git a/MATLAB/Utilities/addCTnoise.m b/MATLAB/Utilities/addCTnoise.m
index abc69131..6a6afc9f 100644
--- a/MATLAB/Utilities/addCTnoise.m
+++ b/MATLAB/Utilities/addCTnoise.m
@@ -1,4 +1,4 @@
-function proj=addCTnoise(proj,varargin)
+function proj=addCTnoise(proj, varargin)
 %ADDCTNOISE adds realistic noise to CT projections
 %   addCTnoise(PROJ):  adds Poisson and Gaussian noise to the input data.
 %
@@ -33,8 +33,8 @@
 %% parse inputs
 
 
-opts=     {'Poisson','Gaussian'};
-defaults= [   1  ,  1 ];
+opts=     {'Poisson','Gaussian', 'Implementation'};
+defaults= [   1  ,  1  ,  1 ];
 
 % Check inputs
 nVarargs = length(varargin);
@@ -85,6 +85,14 @@
                 m=val(1);
                 sigma=val(2);
             end
+        
+        case 'Implementation'
+            if default
+                implementation='matlab';
+            else
+                implementation=val;
+            end
+            
         otherwise
             error('CBCT:addnoise:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_SART_CBCT()']);
     end
@@ -95,20 +103,24 @@
 Im=I0*exp(-proj/max(proj(:)));
 
 % Photon noise + electronic noise
-if areTheseToolboxesInstalled({'MATLAB','Image Processing Toolbox'})
-    Im=imnoise(Im/I0,'poisson')*I0;
-    Im=imnoise(Im/I0,'gaussian',m,sigma/I0)*I0;
-elseif areTheseToolboxesInstalled({'MATLAB','Statistics Toolbox'}) || areTheseToolboxesInstalled({'MATLAB','Statistics and Machine Learning Toolbox'})
-    Im=poissrnd(Im)+randn(size(Im)).*sigma + m;
-else
-    warning(['You don''t have Statistic toolbox, so poisson random noise is not available in MATLAB.',...
-        java.lang.System.getProperty('line.separator').char,...
-        'If you want to add that noise, use the following command:',...
-        java.lang.System.getProperty('line.separator').char,...
-        'Im=poissonrandom(I0*exp(-proj)); Im(Im<0)=1e-6; proj=single(log(I0./Im));',...
-        java.lang.System.getProperty('line.separator').char,...,
-        'With I0 ~ 10000'])
-    Im=randn(size(Im)).*sigma + m; % this one is slower
+if strcmp(implementation, 'matlab')
+    if areTheseToolboxesInstalled({'MATLAB','Image Processing Toolbox'})
+        %disp('Using Image Processing Toolbox');
+        Im=imnoise(Im/I0,'poisson')*I0;
+        Im=imnoise(Im/I0,'gaussian',m,sigma/I0)*I0;
+    elseif areTheseToolboxesInstalled({'MATLAB','Statistics Toolbox'}) || areTheseToolboxesInstalled({'MATLAB','Statistics and Machine Learning Toolbox'})
+        %disp('Using Statistics Toolbox');
+        Im=poissrnd(Im)+randn(size(Im)).*sigma + m;
+    else
+        warning(['You don''t have Image Processing Toolbox nor Statistic Toolbox.',...
+            java.lang.System.getProperty('line.separator').char,...
+            'CUDA version is used.']);
+        %disp('Using CUDA ..');
+        Im=AddNoise(Im, m, sigma);
+    end
+else 
+    %disp('Using CUDA');
+    Im=AddNoise(Im, m, sigma);
 end
 Im(Im<=0)=1e-6;
 proj=single(log(I0./Im))*max(proj(:));
diff --git a/MATLAB/Utilities/cuda_interface/AddNoise.cpp b/MATLAB/Utilities/cuda_interface/AddNoise.cpp
new file mode 100644
index 00000000..3ec22838
--- /dev/null
+++ b/MATLAB/Utilities/cuda_interface/AddNoise.cpp
@@ -0,0 +1,126 @@
+/*-------------------------------------------------------------------------
+ *
+ * MATLAB MEX  functions for Random Number Generator. Check inputs and parses 
+ * MATLAB data to C++ data.
+ *
+ *
+ * CODE by       Ander Biguri, Tomoyuki SADAKANE
+ *
+---------------------------------------------------------------------------
+---------------------------------------------------------------------------
+Copyright (c) 2015, University of Bath and CERN- European Organization for 
+Nuclear Research
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without 
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, 
+this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, 
+this list of conditions and the following disclaimer in the documentation 
+and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors
+may be used to endorse or promote products derived from this software without
+specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
+ ---------------------------------------------------------------------------
+
+Contact: tigre.toolbox@gmail.com
+Codes  : https://github.com/CERN/TIGRE
+--------------------------------------------------------------------------- 
+ */
+
+#include <math.h>
+#include <string.h>
+#include <tmwtypes.h>
+#include <mex.h>
+#include <matrix.h>
+#include <CUDA/RandomNumberGenerator.hpp>
+#include <CUDA/GpuIds.hpp>
+#include <CUDA/gpuUtils.hpp>
+/**
+ * MEX gateway
+ * AddNoise(Im, mu, sigma, "gpuids", gpuids);
+ *   poissrnd(Im)+randn(size(Im)).*sigma + mu;
+ */
+
+void mexFunction(int nlhs, mxArray *plhs[],
+                 int nrhs, mxArray const *prhs[])
+{
+    size_t uiLen = 0;
+    float fGaussMu = 0;
+    float fGaussSigma = 0;
+
+    GpuIds gpuids;
+    if (nrhs==5) {
+        size_t iM = mxGetM(prhs[4]);
+        if (iM != 1) {
+            mexErrMsgIdAndTxt( "CBCT:MEX:RNG:unknown","5th parameter must be a row vector.");
+            return;
+        }
+        size_t uiGpuCount = mxGetN(prhs[4]);
+        if (uiGpuCount == 0) {
+            mexErrMsgIdAndTxt( "CBCT:MEX:RNG:unknown","5th parameter must be a row vector.");
+            return;
+        }
+        int* piGpuIds = (int*)mxGetData(prhs[4]);
+        gpuids.SetIds(uiGpuCount, piGpuIds);
+    } else {
+        int iGpuCount = GetGpuCount();
+        int* piDev = (int*)malloc(iGpuCount * sizeof(int));
+        for (int iI = 0; iI < iGpuCount; ++iI) {
+            piDev[iI] = iI;
+        }
+        gpuids.SetIds(iGpuCount, piDev);
+        free(piDev); piDev = 0;
+    }
+    if (nrhs < 3) {
+        mexErrMsgIdAndTxt("CBCT:CUDA:RNG", "At least three input argumet required.");
+    } else if (nrhs==3 || nrhs==5){
+        size_t mrows = mxGetM(prhs[1]);
+        size_t ncols = mxGetN(prhs[1]);
+        if (mrows!=1 || ncols !=1) {
+            mexErrMsgIdAndTxt("CBCT:CUDA:RNG", "2nd parameter should be 1x1");
+        }
+        mrows = mxGetM(prhs[2]);
+        ncols = mxGetN(prhs[2]);
+        if (mrows!=1 || ncols !=1) {
+            mexErrMsgIdAndTxt("CBCT:CUDA:RNG", "3rd parameter should be 1x1");
+        }
+        fGaussMu    = (float)mxGetScalar(prhs[1]);
+        fGaussSigma = (float)mxGetScalar(prhs[2]);
+    } else if (nrhs>4) {
+        mexErrMsgIdAndTxt("CBCT:CUDA:RNG", "Too many imput argumets");
+    }
+    /////////////// First input argumet.
+    // First input should be an array, whose elements are lambda.
+    mxArray const * const image = prhs[0];
+    float* pfLambdas = static_cast<float*>(mxGetData(image));
+    mwSize const numDims = mxGetNumberOfDimensions(image);  // get dim of image
+    const mwSize *size_img= mxGetDimensions(image); //get size of image
+    uiLen = size_img[0];    // calculate the total length
+    for (int iI = 1; iI < numDims; ++iI) {
+        uiLen *= size_img[iI];
+    }
+    //////////////
+    //prepare outputs
+    // Allocte output image
+    plhs[0] = mxCreateNumericArray(numDims, size_img, mxSINGLE_CLASS, mxREAL);
+    float *imgout =(float*) mxGetPr(plhs[0]);
+    // call CUDA rng
+    poisson_gaussian_1d(pfLambdas, uiLen, fGaussMu, fGaussSigma, imgout, gpuids);
+}

From ed2a85e8cb3abf11733a97bf55d3287fc611f7a8 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Sun, 22 Aug 2021 00:01:52 +0900
Subject: [PATCH 3/7] Add CUDA implementation of random number generator

* Modified and added Python files
---
 Python/setup.py                               | 24 ++++-
 Python/tigre/utilities/CTnoise.py             |  5 +-
 .../cuda_interface/_randomNumberGenerator.pyx | 92 +++++++++++++++++++
 3 files changed, 117 insertions(+), 4 deletions(-)
 create mode 100644 Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx

diff --git a/Python/setup.py b/Python/setup.py
index 1af6d5e4..ce282608 100644
--- a/Python/setup.py
+++ b/Python/setup.py
@@ -450,6 +450,28 @@ def include_headers(filename_list, sdist=False):
     include_dirs=[NUMPY_INCLUDE, CUDA["include"], "../Common/CUDA/"],
 )
 
+
+RandomNumberGenerator_ext = Extension(
+    "_RandomNumberGenerator",
+    sources=include_headers(
+        [
+            "../Common/CUDA/TIGRE_common.cpp",
+            "../Common/CUDA/RandomNumberGenerator.cu",
+            "../Common/CUDA/GpuIds.cpp",
+            "../Common/CUDA/gpuUtils.cu",
+            "tigre/utilities/cuda_interface/_randomNumberGenerator.pyx",
+        ],
+        sdist=sys.argv[1] == "sdist",
+    ),
+    define_macros=[("IS_FOR_PYTIGRE", None)],
+    library_dirs=[CUDA["lib64"]],
+    libraries=["cudart"],
+    language="c++",
+    runtime_library_dirs=[CUDA["lib64"]] if not IS_WINDOWS else None,
+    include_dirs=[NUMPY_INCLUDE, CUDA["include"], "../Common/CUDA/"],
+)
+
+
 setup(
     name="pytigre",
     version="2.1.0",
@@ -457,7 +479,7 @@ def include_headers(filename_list, sdist=False):
     packages=find_packages(),
     include_package_data=True,
     data_files=[("data", ["../Common/data/head.mat"])],
-    ext_modules=[Ax_ext, Atb_ext, tvdenoising_ext, minTV_ext, AwminTV_ext, gpuUtils_ext],
+    ext_modules=[Ax_ext, Atb_ext, tvdenoising_ext, minTV_ext, AwminTV_ext, gpuUtils_ext, RandomNumberGenerator_ext],
     py_modules=["tigre.py"],
     cmdclass={"build_ext": BuildExtension},
     install_requires=["Cython", "matplotlib", "numpy", "scipy"],
diff --git a/Python/tigre/utilities/CTnoise.py b/Python/tigre/utilities/CTnoise.py
index 2c13ac83..0706e7c5 100644
--- a/Python/tigre/utilities/CTnoise.py
+++ b/Python/tigre/utilities/CTnoise.py
@@ -1,5 +1,5 @@
 import numpy as np
-
+import _RandomNumberGenerator as RNG
 
 def add(projections, Gaussian=None, Poisson=None):
 
@@ -21,8 +21,7 @@ def add(projections, Gaussian=None, Poisson=None):
     max_proj = np.max(projections)
     projections = Poisson * np.exp(-projections / max_proj)
 
-    projections = np.random.poisson(projections)
-    projections = projections + np.random.normal(Gaussian[0], Gaussian[1], size=projections.shape)
+    projections = RNG.add_noise(projections, Gaussian[0], Gaussian[1])
 
     projections = -np.log(projections / Poisson) * max_proj
     projections = np.float32(projections)
diff --git a/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx b/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx
new file mode 100644
index 00000000..9f684b6e
--- /dev/null
+++ b/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx
@@ -0,0 +1,92 @@
+#   This file is part of the TIGRE Toolbox
+
+#   Copyright (c) 2015, University of Bath and
+#                       CERN-European Organization for Nuclear Research
+#                       All rights reserved.
+
+#   License:            Open Source under BSD.
+#                       See the full license at
+#                       https://github.com/CERN/TIGRE/license.txt
+
+#   Contact:            tigre.toolbox@gmail.com
+#   Codes:              https://github.com/CERN/TIGRE/
+# --------------------------------------------------------------------------
+#   Coded by:          Ander Biguri, Tomoyuki SADAKANE
+
+cimport numpy as np 
+import numpy as np
+from tigre.utilities.errors import TigreCudaCallError
+from tigre.utilities.cuda_interface._gpuUtils cimport GpuIds as c_GpuIds, convert_to_c_gpuids, free_c_gpuids
+
+np.import_array()
+
+from libc.stdlib cimport malloc, free 
+
+cdef extern from "numpy/arrayobject.h":
+    void PyArray_ENABLEFLAGS(np.ndarray arr, int flags)
+    void PyArray_CLEARFLAGS(np.ndarray arr, int flags)
+
+cdef extern from "RandomNumberGenerator.hpp":
+    cdef void poisson_1d(float* img, size_t uiLen, float* dst, c_GpuIds gpuids)
+    cdef void poisson_gaussian_1d(float* img, size_t uiLen, float fGaussMu, float fGaussSigma, float* dst, c_GpuIds gpuids)
+
+def cuda_raise_errors(error_code):
+    if error_code:
+        raise TigreCudaCallError('RandomNumberGenerator::', error_code)
+
+def add_poisson(np.ndarray[np.float32_t, ndim=3] src, gpuids=None):
+    # print("add_poisson()")
+    cdef c_GpuIds* c_gpuids = convert_to_c_gpuids(gpuids)
+    if not c_gpuids:
+        raise MemoryError()
+
+    cdef np.npy_intp size_img[3]
+    size_img[0]= <np.npy_intp> src.shape[0]
+    size_img[1]= <np.npy_intp> src.shape[1]
+    size_img[2]= <np.npy_intp> src.shape[2]
+
+    cdef float* c_imgout = <float*> malloc(size_img[0] *size_img[1] *size_img[2]* sizeof(float))
+
+    cdef long imgsize[3]
+    imgsize[0] = <long> size_img[2]
+    imgsize[1] = <long> size_img[1]
+    imgsize[2] = <long> size_img[0]
+
+    cdef float* c_src = <float*> src.data
+    cdef np.npy_intp c_uiSigLen = <np.npy_intp> (size_img[0] *size_img[1] *size_img[2])
+    cuda_raise_errors(poisson_1d(c_src, c_uiSigLen, c_imgout, c_gpuids[0]))
+    imgout = np.PyArray_SimpleNewFromData(3, size_img, np.NPY_FLOAT32, c_imgout)
+    PyArray_ENABLEFLAGS(imgout, np.NPY_OWNDATA)
+    
+    return imgout
+
+def add_noise(np.ndarray[np.float32_t, ndim=3] poisson_lambda,
+              np.float32_t gaussian_mu,
+              np.float32_t gaussian_sigma,
+              gpuids=None):
+    # print("add_noise()")
+    cdef c_GpuIds* c_gpuids = convert_to_c_gpuids(gpuids)
+    if not c_gpuids:
+        raise MemoryError()
+
+    cdef np.npy_intp size_img[3]
+    size_img[0]= <np.npy_intp> poisson_lambda.shape[0]
+    size_img[1]= <np.npy_intp> poisson_lambda.shape[1]
+    size_img[2]= <np.npy_intp> poisson_lambda.shape[2]
+
+    cdef float* c_imgout = <float*> malloc(size_img[0] *size_img[1] *size_img[2]* sizeof(float))
+
+    cdef long imgsize[3]
+    imgsize[0] = <long> size_img[2]
+    imgsize[1] = <long> size_img[1]
+    imgsize[2] = <long> size_img[0]
+
+    cdef float* c_src = <float*> poisson_lambda.data
+    cdef np.npy_intp c_uiSigLen = <np.npy_intp> (size_img[0] *size_img[1] *size_img[2])
+    cdef np.float32_t c_fGaussMu = gaussian_mu
+    cdef np.float32_t c_fGaussSigma = gaussian_sigma
+    cuda_raise_errors(poisson_gaussian_1d(c_src, c_uiSigLen, c_fGaussMu, c_fGaussSigma, c_imgout, c_gpuids[0]))
+    imgout = np.PyArray_SimpleNewFromData(3, size_img, np.NPY_FLOAT32, c_imgout)
+    PyArray_ENABLEFLAGS(imgout, np.NPY_OWNDATA)
+    
+    return imgout

From c8c4b676b2fbe26849285f86561420dd0a9f8546 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Tue, 31 Aug 2021 21:43:15 +0900
Subject: [PATCH 4/7] Make the default of MATLAB addCTnoise CUDA

* The default implementation is changed to CUDA.
* See #320, #321
---
 MATLAB/Utilities/addCTnoise.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/MATLAB/Utilities/addCTnoise.m b/MATLAB/Utilities/addCTnoise.m
index 6a6afc9f..3c44cb35 100644
--- a/MATLAB/Utilities/addCTnoise.m
+++ b/MATLAB/Utilities/addCTnoise.m
@@ -7,7 +7,7 @@
 %   'Poisson' : changes the average photon count for Poisson noise. Default
 %               value is 1e5
 %   'Gaussian': changes the mean and standard deviation of the electronic
-%               Gaussian noise. Default value is [0 10]
+%               Gaussian noise. Default value is [0 0.5]
 %
 % The computation of the noise projection follows
 % Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction
@@ -88,7 +88,7 @@
         
         case 'Implementation'
             if default
-                implementation='matlab';
+                implementation='cuda';
             else
                 implementation=val;
             end
@@ -118,7 +118,7 @@
         %disp('Using CUDA ..');
         Im=AddNoise(Im, m, sigma);
     end
-else 
+else % if implementation == 'cuda'
     %disp('Using CUDA');
     Im=AddNoise(Im, m, sigma);
 end

From 0e89d5cda745e1fe18d7ba5f86d3baf6ba67dba3 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Thu, 2 Sep 2021 22:43:36 +0900
Subject: [PATCH 5/7] Update the author name

---
 Common/CUDA/RandomNumberGenerator.cu                           | 3 +--
 Common/CUDA/RandomNumberGenerator.hpp                          | 3 +--
 MATLAB/Utilities/cuda_interface/AddNoise.cpp                   | 2 +-
 .../tigre/utilities/cuda_interface/_randomNumberGenerator.pyx  | 2 +-
 4 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/Common/CUDA/RandomNumberGenerator.cu b/Common/CUDA/RandomNumberGenerator.cu
index f9cd2188..d7d1224a 100644
--- a/Common/CUDA/RandomNumberGenerator.cu
+++ b/Common/CUDA/RandomNumberGenerator.cu
@@ -4,8 +4,7 @@
  *
  * Adds noise of Poisson and normal distribution to the input.
  *
- * CODE by       Ander Biguri
- *               Tomoyuki SADAKANE
+ * CODE by       Tomoyuki SADAKANE
  * ---------------------------------------------------------------------------
  * ---------------------------------------------------------------------------
  * Copyright (c) 2015, University of Bath and CERN- European Organization for
diff --git a/Common/CUDA/RandomNumberGenerator.hpp b/Common/CUDA/RandomNumberGenerator.hpp
index d6f6c6c2..4ba68d8d 100644
--- a/Common/CUDA/RandomNumberGenerator.hpp
+++ b/Common/CUDA/RandomNumberGenerator.hpp
@@ -4,8 +4,7 @@
  *
  * Adds noise of Poisson and normal distribution to the input.
  *
- * CODE by       Ander Biguri
- *               Tomoyuki SADAKANE
+ * CODE by       Tomoyuki SADAKANE
  * ---------------------------------------------------------------------------
  * ---------------------------------------------------------------------------
  * Copyright (c) 2015, University of Bath and CERN- European Organization for
diff --git a/MATLAB/Utilities/cuda_interface/AddNoise.cpp b/MATLAB/Utilities/cuda_interface/AddNoise.cpp
index 3ec22838..e12d11e5 100644
--- a/MATLAB/Utilities/cuda_interface/AddNoise.cpp
+++ b/MATLAB/Utilities/cuda_interface/AddNoise.cpp
@@ -4,7 +4,7 @@
  * MATLAB data to C++ data.
  *
  *
- * CODE by       Ander Biguri, Tomoyuki SADAKANE
+ * CODE by       Tomoyuki SADAKANE
  *
 ---------------------------------------------------------------------------
 ---------------------------------------------------------------------------
diff --git a/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx b/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx
index 9f684b6e..0e4da405 100644
--- a/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx
+++ b/Python/tigre/utilities/cuda_interface/_randomNumberGenerator.pyx
@@ -11,7 +11,7 @@
 #   Contact:            tigre.toolbox@gmail.com
 #   Codes:              https://github.com/CERN/TIGRE/
 # --------------------------------------------------------------------------
-#   Coded by:          Ander Biguri, Tomoyuki SADAKANE
+#   Coded by:          Tomoyuki SADAKANE
 
 cimport numpy as np 
 import numpy as np

From 8ca3d19aaacca632a000266a2e280f4c5d680219 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Thu, 2 Sep 2021 22:45:19 +0900
Subject: [PATCH 6/7] Fix typos

---
 MATLAB/Utilities/addCTnoise.m | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/MATLAB/Utilities/addCTnoise.m b/MATLAB/Utilities/addCTnoise.m
index 3c44cb35..cfd55038 100644
--- a/MATLAB/Utilities/addCTnoise.m
+++ b/MATLAB/Utilities/addCTnoise.m
@@ -71,7 +71,7 @@
                     I0=max(proj(:))/5;
                 end
             else
-                if ~isscalar(val);error('CBCT:addnoise:WorngInput','Input to Poisson should be scalar');end
+                if ~isscalar(val);error('CBCT:addnoise:WrongInput','Input to Poisson should be scalar');end
                 I0=val;
             end
             
@@ -81,7 +81,7 @@
                 m=0;
                 sigma=0.5;
             else
-                if (size(val,2)~=2);error('CBCT:addnoise:WorngInput','Input to Gaussian should be 1x2');end
+                if (size(val,2)~=2);error('CBCT:addnoise:WrongInput','Input to Gaussian should be 1x2');end
                 m=val(1);
                 sigma=val(2);
             end
@@ -94,7 +94,7 @@
             end
             
         otherwise
-            error('CBCT:addnoise:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in OS_SART_CBCT()']);
+            error('CBCT:addnoise:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in addCTnoise()']);
     end
 end
 

From 8dcca6db517cfcd4d336a8e5f0b0c26e6b033518 Mon Sep 17 00:00:00 2001
From: Tomoyuki SADAKANE <40597344+tsadakane@users.noreply.github.com>
Date: Thu, 2 Sep 2021 22:47:34 +0900
Subject: [PATCH 7/7] Add to Changelog.md the feature description

---
 changelog.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/changelog.md b/changelog.md
index ae430115..f4bd95c4 100644
--- a/changelog.md
+++ b/changelog.md
@@ -14,6 +14,7 @@ Many contributors have made this possible. Sincere thanks.
 - Code quality optional code added, to be able to maintain decent python.
 - Demos for python available
 - Add 3D Shepp-Logan head phantom to Python.
+- Replace noise generator with CUDA implementation.
 
 ## BugFix