diff --git a/Common/CUDA/PICCS.cu b/Common/CUDA/PICCS.cu
index 481ede08..0ddae170 100644
--- a/Common/CUDA/PICCS.cu
+++ b/Common/CUDA/PICCS.cu
@@ -54,7 +54,7 @@ Codes  : https://github.com/CERN/TIGRE
 #define MAXTHREADS 1024
 
 #include "PICCS.hpp"
-
+#include "gpuUtils.hpp"
 
 
 
@@ -114,17 +114,21 @@ do { \
         unsigned long size2d = rows*cols;
         unsigned long long idx = z * size2d + y * cols + x;
         
-        float uidx = u[idx];
+        // float uidx = u[idx];
+        float uidx = 0;
         
         if ( z - 1 >= 0 && z<depth) {
+            uidx = u[idx];
             grad[0] = (uidx-u[(z-1)*size2d + y*cols + x]) ;
         }
-        
-        if ( y - 1 >= 0 && y<rows){
+        //z may be out of bounds, need to add a condition check z < depth.
+        if ( y - 1 >= 0 && y<rows && z<depth){
+            uidx = u[idx];
             grad[1] = (uidx-u[z*size2d + (y-1)*cols + x]) ;
         }
-        
-        if ( x - 1 >= 0 && x<cols) {
+        //z may be out of bounds, need to add a condition check z < depth.
+        if ( x - 1 >= 0 && x<cols && z<depth) {
+            uidx = u[idx];
             grad[2] = (uidx-u[z*size2d + y*cols + (x-1)]);
         }
     }
diff --git a/Common/CUDA/PICCS.hpp b/Common/CUDA/PICCS.hpp
index e3592dbb..7f7d6763 100644
--- a/Common/CUDA/PICCS.hpp
+++ b/Common/CUDA/PICCS.hpp
@@ -55,7 +55,7 @@ Codes  : https://github.com/CERN/TIGRE
 #include "TIGRE_common.hpp"
 #include "GpuIds.hpp"
 
-void piccs_tv(const float* img,const float* prior, float* dst,float alpha, float ratio, const long* image_size, int maxIter, const GpuIds& gpuids);
+void piccs_tv(float* img,float* prior, float* dst,float alpha, float ratio, const long* image_size, int maxIter, const GpuIds& gpuids);
 
 
 #endif
\ No newline at end of file
diff --git a/Python/setup.py b/Python/setup.py
index 40bc3b3f..60586a3f 100644
--- a/Python/setup.py
+++ b/Python/setup.py
@@ -464,6 +464,25 @@ def include_headers(filename_list, sdist=False):
     include_dirs=[NUMPY_INCLUDE, CUDA["include"], "../Common/CUDA/"],
 )
 
+PICCS_ext = Extension(
+    "_PICCS",
+    sources=include_headers(
+        [
+            "../Common/CUDA/TIGRE_common.cpp",
+            "../Common/CUDA/PICCS.cu",
+            "../Common/CUDA/GpuIds.cpp",
+            "../Common/CUDA/gpuUtils.cu",
+            "tigre/utilities/cuda_interface/_PICCS.pyx",
+        ],
+        sdist=sys.argv[1] == "sdist",
+    ),
+    define_macros=define_macros,
+    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/"],
+)
 
 gpuUtils_ext = Extension(
     "_gpuUtils",
@@ -510,7 +529,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, tv_proximal_ext, minTV_ext, AwminTV_ext, gpuUtils_ext, RandomNumberGenerator_ext],
+    ext_modules=[Ax_ext, Atb_ext, tv_proximal_ext, minTV_ext, AwminTV_ext, PICCS_ext, gpuUtils_ext, RandomNumberGenerator_ext],
     py_modules=["tigre.py"],
     cmdclass={"build_ext": BuildExtension},
     install_requires=["Cython", "matplotlib", "numpy", "scipy", "tqdm"],
diff --git a/Python/tigre/algorithms/iterative_recon_alg.py b/Python/tigre/algorithms/iterative_recon_alg.py
index fde3bfef..560ce3db 100644
--- a/Python/tigre/algorithms/iterative_recon_alg.py
+++ b/Python/tigre/algorithms/iterative_recon_alg.py
@@ -6,6 +6,7 @@
 import numpy as np
 from _AwminTV import AwminTV
 from _minTV import minTV
+from _PICCS import PICCS
 from tigre.algorithms.single_pass_algorithms import FDK
 from tigre.utilities.Atb import Atb
 from tigre.utilities.Ax import Ax
@@ -355,6 +356,11 @@ def minimizeAwTV(self, res_prev, dtvg):
             self.gpuids = GpuIds()
         return AwminTV(res_prev, dtvg, self.numiter_tv, self.delta, self.gpuids)
 
+    def PICCS(self, res_prev, res_prior, dtvg, ratio):
+        if self.gpuids is None:
+            self.gpuids = GpuIds()
+        return PICCS(res_prev, res_prior, dtvg, ratio, self.numiter_tv, self.gpuids)
+    
     def error_measurement(self, res_prev, iter):
         if self.Quameasopts is not None:
             self.lq[:, iter] = MQ(self.res, res_prev, self.Quameasopts)
diff --git a/Python/tigre/algorithms/pocs_algorithms.py b/Python/tigre/algorithms/pocs_algorithms.py
index 2ea51f1f..52ddb922 100644
--- a/Python/tigre/algorithms/pocs_algorithms.py
+++ b/Python/tigre/algorithms/pocs_algorithms.py
@@ -112,6 +112,8 @@ def __init__(self, proj, geo, angles, niter, **kwargs):
         self.alpha = 0.002 if "alpha" not in kwargs else kwargs["alpha"]
         self.alpha_red = 0.95 if "alpha_red" not in kwargs else kwargs["alpha_red"]
         self.rmax = 0.95 if "rmax" not in kwargs else kwargs["rmax"]
+        self.ratio = 0.5 if "ratio" not in kwargs else kwargs["ratio"]
+        self.regularization = "minimizeTV" if "regularization" not in kwargs else kwargs["regularization"]
         if "maxl2err" not in kwargs:
             self.epsilon = (
                 im3DNORM(Ax(FDK(proj, geo, angles, gpuids=self.gpuids), geo, angles) - proj, 2)
@@ -128,6 +130,10 @@ def __init__(self, proj, geo, angles, niter, **kwargs):
     def run_main_iter(self):
         stop_criteria = False
         n_iter = 0
+        if (self.regularization=="PICCS"):
+            # Modify the `set_res` function in the `Python\tigre\algorithms\iterative_recon_alg.py` file as needed.
+            res_prior = copy.deepcopy(self.res)
+
         while not stop_criteria:
             if self.verbose:
                 self._estimate_time_until_completion(n_iter)
@@ -146,7 +152,12 @@ def run_main_iter(self):
                 dtvg = self.alpha * dp
 
             res_prev = copy.deepcopy(self.res)
-            self.res = getattr(self, self.regularization)(self.res, dtvg)
+            # self.res = getattr(self, self.regularization)(self.res, dtvg)
+            if (self.regularization=="PICCS"):
+                self.res = getattr(self, self.regularization)(self.res, res_prior, dtvg, self.ratio)
+            else:
+                self.res = getattr(self, self.regularization)(self.res, dtvg)
+            
             dg_vec = self.res - res_prev
             dg = im3DNORM(dg_vec, 2)
 
diff --git a/Python/tigre/utilities/cuda_interface/_PICCS.pyx b/Python/tigre/utilities/cuda_interface/_PICCS.pyx
new file mode 100644
index 00000000..aab5f6c9
--- /dev/null
+++ b/Python/tigre/utilities/cuda_interface/_PICCS.pyx
@@ -0,0 +1,50 @@
+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 "PICCS.hpp":
+    cdef void piccs_tv(float* img, float* prior, float* dst, float alpha, float ratio, long* image_size, int maxiter, c_GpuIds gpuids)
+
+
+def cuda_raise_errors(error_code):
+    if error_code:
+        raise TigreCudaCallError('PICCS:PICCS_TV:', error_code)
+
+
+def PICCS(np.ndarray[np.float32_t, ndim=3] src,np.ndarray[np.float32_t, ndim=3] prior,float alpha = 15.0,float ratio=0.5,int maxiter = 100, gpuids=None):
+    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(<unsigned long>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]
+    
+    src = np.ascontiguousarray(src)
+    prior = np.ascontiguousarray(prior)
+
+    cdef float* c_src = <float*> src.data
+    cdef float* c_prior = <float*> prior.data
+    cdef np.npy_intp c_maxiter = <np.npy_intp> maxiter
+    cuda_raise_errors(piccs_tv(c_src, c_prior, c_imgout, alpha, ratio, imgsize, c_maxiter, c_gpuids[0]))
+    imgout = np.PyArray_SimpleNewFromData(3, size_img, np.NPY_FLOAT32, c_imgout)
+    PyArray_ENABLEFLAGS(imgout, np.NPY_OWNDATA)
+    
+    return imgout