Skip to content

Modify and update the code related to the PICCS algorithm module. #650

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Modify and update the code related to the PICCS algorithm module. #650

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dcf1883
Modify and update the code related to the PICCS algorithm module.
Apr 1, 2025
7a37a9e
little change in the PICCS module
Apr 1, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 10 additions & 6 deletions Common/CUDA/PICCS.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ Codes : https://github.com/CERN/TIGRE
#define MAXTHREADS 1024

#include "PICCS.hpp"

#include "gpuUtils.hpp"



Expand Down Expand Up @@ -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)]);
}
}
Expand Down
2 changes: 1 addition & 1 deletion Common/CUDA/PICCS.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did you remove const here?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This modification is an attempt I made while debugging the incorrect return value of this operator. Keeping const is a better choice.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@coldly01 makes sense :) can you put it back then?



#endif
21 changes: 20 additions & 1 deletion Python/setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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"],
Expand Down
6 changes: 6 additions & 0 deletions Python/tigre/algorithms/iterative_recon_alg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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_TV(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)
Expand Down
13 changes: 12 additions & 1 deletion Python/tigre/algorithms/pocs_algorithms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hum not sure about this. res is the current variable right? You may want to initialize the algorithm to one image, but have a prior of another, so I think there should be a new kwarg for prior

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My approach is to add the following code in the set_res function.

if init == "PICCS_priorimage":
                file_path='./fdk/fdk_360ang_test.raw'
                data = np.fromfile(file_path, dtype=np.float32)
                nVoxel = self.geo.nVoxel
                self.res=data.reshape(nVoxel)

And then call os_asd_pocs in demo.py like this.

os_asd_pocs_output=algs.os_asd_pocs(proj, geo, angles, niter, regularization="PICCS", alpha=piccs_alpha, ratio=piccs_ratio, init="PICCS_priorimage", gpuids=gpuids)

Of course a better approach would be to add a kwarg with a name like prior_image_filepath

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that should be the right way.

In MATLAB we have PICCS as a separate function, ideally we want to have an alias (like os_asd_pocs) for piccs and os-piccs that force you to have a kwarg of prior_image, but not a filepath, but a volume of the same shape as geo.nVoxel.

Let me know if you want to implement this change, or want me to do it before the merge (I will make commits in your branch if so).


while not stop_criteria:
if self.verbose:
self._estimate_time_until_completion(n_iter)
Expand All @@ -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)

Expand Down
50 changes: 50 additions & 0 deletions Python/tigre/utilities/cuda_interface/_PICCS.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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