Skip to content

Commit f4a12d9

Browse files
dkazancdkazanc
authored
Merge pull request #252 from DiamondLightSource/cupy-update-radway-65
Pyproject updates for Cupy 13-14
2 parents ea0718f + 660ec2d commit f4a12d9

7 files changed

Lines changed: 63 additions & 172 deletions

File tree

‎README.rst‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ Some of the methods also have been optimised to ensure higher computational effi
99
The purpose of HTTomolibGPU
1010
===========================
1111

12-
Although **HTTomolibGPU** can be used as a stand-alone library, it has been specifically developed to work together with the
12+
Although **HTTomolibGPU** can be used as a stand-alone library, it has been specifically developed to work together with the
1313
`HTTomo <https://diamondlightsource.github.io/httomo/>`_ package as
1414
its backend for data processing. HTTomo is a user interface (UI) written in Python for fast big tomographic data processing using
1515
MPI protocols or as well serially.
@@ -34,7 +34,7 @@ Conda environment
3434
3535
$ conda create --name httomolibgpu # create a fresh conda environment
3636
$ conda activate httomolibgpu # activate the environment
37-
$ conda install conda-forge::cupy==12.3.0
37+
$ conda install conda-forge::cupy==14.0.*
3838
$ pip install httomolibgpu
3939
4040
Setup the development environment:
@@ -43,6 +43,6 @@ Setup the development environment:
4343
.. code-block:: console
4444
4545
$ git clone git@github.com:DiamondLightSource/httomolibgpu.git # clone the repo
46-
$ conda env create --name httomolibgpu -c conda-forge cupy==12.3.0 # install dependencies
46+
$ conda env create --name httomolibgpu -c conda-forge cupy==14.0.* # install dependencies
4747
$ conda activate httomolibgpu # activate the environment
4848
$ pip install -e ./httomolibgpu[dev] # editable/development mode

‎httomolibgpu/cuda_kernels/__init__.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,5 +25,5 @@ def load_cuda_module(
2525
code += f.read()
2626

2727
return cp.RawModule(
28-
options=("-std=c++11", *options), code=code, name_expressions=name_expressions
28+
options=("-std=c++17", *options), code=code, name_expressions=name_expressions
2929
)

‎httomolibgpu/prep/normalize.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ def dark_flat_field_correction(
123123
"float32 out",
124124
kernel,
125125
kernel_name,
126-
options=("-std=c++11",),
126+
options=("-std=c++17",),
127127
loop_prep="constexpr float eps = 1.0e-07;",
128128
no_return=True,
129129
)

‎httomolibgpu/prep/phase.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -207,7 +207,7 @@ def paganin_filter(
207207
mem_stack.free(ifft_plan.work_area.mem.size)
208208
else:
209209
with ifft_plan:
210-
ifft_filtered_data = ifft2(fft_data, axes=(-2, -1), overwrite_x=True).real
210+
ifft_filtered_data = ifft2(fft_data, axes=(-2, -1), overwrite_x=True)
211211
del fft_data
212212
del ifft_plan
213213
del ifft_input
@@ -230,6 +230,7 @@ def paganin_filter(
230230
return mem_stack.highwater
231231

232232
# crop the padded filtered data:
233+
# .real should not be used for now, see: https://github.com/cupy/cupy/issues/9750
233234
data = ifft_filtered_data[slc_indices].astype(cp.float32)
234235
del ifft_filtered_data
235236

‎httomolibgpu/recon/algorithm.py‎

Lines changed: 49 additions & 159 deletions
Original file line numberDiff line numberDiff line change
@@ -405,12 +405,12 @@ def SIRT3d_tomobar(
405405
center,
406406
detector_pad,
407407
recon_size,
408+
1,
408409
gpu_id,
409-
datafidelity="LS",
410410
)
411411

412412
_data_ = {
413-
"projection_norm_data": data,
413+
"projection_data": data,
414414
"data_axes_labels_order": input_data_axis_labels,
415415
} # data dictionary
416416
_algorithm_ = {
@@ -484,11 +484,17 @@ def CGLS3d_tomobar(
484484
###################################
485485

486486
RecToolsCP = _instantiate_iterative_recon_class(
487-
data, angles, center, detector_pad, recon_size, gpu_id, datafidelity="LS"
487+
data,
488+
angles,
489+
center,
490+
detector_pad,
491+
recon_size,
492+
1,
493+
gpu_id,
488494
)
489495

490496
_data_ = {
491-
"projection_norm_data": data,
497+
"projection_data": data,
492498
"data_axes_labels_order": input_data_axis_labels,
493499
} # data dictionary
494500
_algorithm_ = {
@@ -511,6 +517,7 @@ def FISTA3d_tomobar(
511517
recon_mask_radius: Optional[float] = 0.95,
512518
iterations: int = 20,
513519
subsets_number: int = 6,
520+
data_fidelity: str = "LS",
514521
regularisation_type: Literal["ROF_TV", "PD_TV"] = "PD_TV",
515522
regularisation_parameter: float = 0.000001,
516523
regularisation_iterations: int = 50,
@@ -544,6 +551,8 @@ def FISTA3d_tomobar(
544551
The number of FISTA algorithm iterations.
545552
subsets_number: int
546553
The number of the ordered subsets to accelerate convergence. Keep the value bellow 10 to avoid divergence.
554+
data_fidelity: str
555+
Data fidelity given as 'LS' (Least Squares), 'PWLS' (Penalised Weightes LS).
547556
regularisation_type: str
548557
A method to use for regularisation. Currently PD_TV and ROF_TV are available.
549558
regularisation_parameter: float
@@ -586,19 +595,25 @@ def FISTA3d_tomobar(
586595
###################################
587596

588597
RecToolsCP = _instantiate_iterative_recon_class(
589-
data, angles, center, detector_pad, recon_size, gpu_id, datafidelity="LS"
598+
data,
599+
angles,
600+
center,
601+
detector_pad,
602+
recon_size,
603+
subsets_number,
604+
gpu_id,
590605
)
591606

592607
_data_ = {
593-
"projection_norm_data": data,
594-
"OS_number": subsets_number,
608+
"projection_data": data,
609+
"data_fidelity": data_fidelity,
595610
"data_axes_labels_order": input_data_axis_labels,
596611
}
597612
lc = RecToolsCP.powermethod(_data_) # calculate Lipschitz constant (run once)
598613

599614
_algorithm_ = {
600615
"iterations": iterations,
601-
"lipschitz_const": lc.get(),
616+
"lipschitz_const": lc,
602617
"nonnegativity": nonnegativity,
603618
"recon_mask_radius": recon_mask_radius,
604619
}
@@ -625,7 +640,8 @@ def ADMM3d_tomobar(
625640
recon_mask_radius: Optional[float] = 0.95,
626641
iterations: int = 3,
627642
subsets_number: int = 24,
628-
initialisation: Literal["FBP", "CGLS", None] = "FBP",
643+
data_fidelity: str = "LS",
644+
initialisation: Literal["FBP", "CGLS", "SIRT", None] = "FBP",
629645
ADMM_rho_const: float = 1.0,
630646
ADMM_relax_par: float = 1.7,
631647
regularisation_type: str = "PD_TV",
@@ -662,9 +678,11 @@ def ADMM3d_tomobar(
662678
more than 10 without. Assuming that the subsets_number is reasonably large (>12).
663679
subsets_number: int
664680
The number of the ordered subsets to accelerate convergence. The recommended range is between 12 to 24.
681+
data_fidelity: str
682+
Data fidelity given as 'LS' (Least Squares), 'PWLS' (Penalised Weightes LS).
665683
initialisation: str, optional
666-
Initialise ADMM with the reconstructed image to reduce the number of iterations and accelerate. Choose between 'CGLS' when data
667-
is noisy and undersampled, 'FBP' when data is of better quality (default) or None.
684+
Initialise ADMM with the reconstructed image to reduce the number of iterations and accelerate. Choose between 'CGLS' or 'SIRT' when data
685+
is noisy and/or undersampled. Choose 'FBP' when the data is of better quality (default) or None.
668686
ADMM_rho_const: float
669687
Convergence related parameter for ADMM, higher values lead to slower convergence, but too small values can destabilise the iterations.
670688
Recommended range is between 0.9 and 2.0.
@@ -712,7 +730,7 @@ def ADMM3d_tomobar(
712730
initialisation,
713731
[str, type(None)],
714732
"initialisation",
715-
["FBP", "CGLS"],
733+
["FBP", "CGLS", "SIRT"],
716734
methods_name,
717735
)
718736
__check_variable_type(ADMM_rho_const, [float], "ADMM_rho_const", [], methods_name)
@@ -759,152 +777,18 @@ def ADMM3d_tomobar(
759777
),
760778
requirements="C",
761779
)
762-
else:
763-
initialisation_vol = None
764-
765-
RecToolsCP = _instantiate_iterative_recon_class(
766-
data, angles, center, detector_pad, recon_size, gpu_id, datafidelity="LS"
767-
)
768-
769-
_data_ = {
770-
"projection_norm_data": data,
771-
"OS_number": subsets_number,
772-
"data_axes_labels_order": input_data_axis_labels,
773-
}
774-
775-
_algorithm_ = {
776-
"initialise": initialisation_vol,
777-
"iterations": iterations,
778-
"nonnegativity": nonnegativity,
779-
"recon_mask_radius": recon_mask_radius,
780-
"ADMM_rho_const": ADMM_rho_const,
781-
"ADMM_relax_par": ADMM_relax_par,
782-
}
783-
784-
_regularisation_ = {
785-
"method": regularisation_type, # Selected regularisation method
786-
"regul_param": regularisation_parameter, # Regularisation parameter
787-
"iterations": regularisation_iterations, # The number of regularisation iterations
788-
"half_precision": regularisation_half_precision, # enabling half-precision calculation
789-
}
790-
791-
reconstruction = RecToolsCP.ADMM(_data_, _algorithm_, _regularisation_)
792-
cp._default_memory_pool.free_all_blocks()
793-
return cp.require(cp.swapaxes(reconstruction, 0, 1), requirements="C")
794-
795-
796-
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ##
797-
def _instantiate_direct_recon_class(
798-
data: cp.ndarray,
799-
angles: np.ndarray,
800-
center: Optional[float] = None,
801-
detector_pad: Union[bool, int] = False,
802-
recon_size: Optional[int] = None,
803-
recon_mask_radius: float = 0.95,
804-
iterations: int = 3,
805-
subsets_number: int = 24,
806-
initialisation: Optional[str] = "FBP",
807-
ADMM_rho_const: float = 1.0,
808-
ADMM_relax_par: float = 1.7,
809-
regularisation_type: str = "PD_TV",
810-
regularisation_parameter: float = 0.0025,
811-
regularisation_iterations: int = 40,
812-
regularisation_half_precision: bool = True,
813-
nonnegativity: bool = False,
814-
gpu_id: int = 0,
815-
) -> cp.ndarray:
816-
"""
817-
An Alternating Direction Method of Multipliers method with various types of regularisation or
818-
denoising operations :cite:`kazantsev2019ccpi` (currently accepts ROF_TV and PD_TV regularisations only).
819-
For more information see :ref:`_method_ADMM3d_tomobar`.
820-
821-
Parameters
822-
----------
823-
data : cp.ndarray
824-
Projection data as a CuPy array.
825-
angles : np.ndarray
826-
An array of angles given in radians.
827-
center : float, optional
828-
The center of rotation (CoR).
829-
detector_pad : bool, int
830-
Detector width padding with edge values to remove circle/arc type artifacts in the reconstruction. Set to True to perform
831-
an automated padding or specify a certain value as an integer.
832-
recon_size : int, optional
833-
The [recon_size, recon_size] shape of the reconstructed slice in pixels.
834-
By default (None), the reconstructed size will be the dimension of the horizontal detector.
835-
recon_mask_radius: float
836-
The radius of the circular mask that applies to the reconstructed slice in order to crop
837-
out some undesirable artifacts. The values outside the given diameter will be set to zero.
838-
To implement the cropping one can use the range [0.7-1.0] or set to 2.0 when no cropping required.
839-
iterations : int
840-
The number of ADMM algorithm iterations. The recommended range is between 3 to 5 with initialisation and
841-
more than 10 without. Assuming that the subsets_number is reasonably large (>12).
842-
subsets_number: int
843-
The number of the ordered subsets to accelerate convergence. The recommended range is between 12 to 24.
844-
initialisation: str, optional
845-
Initialise ADMM with the reconstructed image to reduce the number of iterations and accelerate. Choose between 'CGLS' when data
846-
is noisy and undersampled, 'FBP' when data is of better quality (default) or None.
847-
ADMM_rho_const: float
848-
Convergence related parameter for ADMM, higher values lead to slower convergence, but too small values can destabilise the iterations.
849-
Recommended range is between 0.9 and 2.0.
850-
ADMM_relax_par: Relaxation parameter which can lead to acceleration of the algorithm, keep it in the range between 1.5 and 1.8 to avoid divergence. regularisation_type: str
851-
A method to use for regularisation. Currently PD_TV and ROF_TV are available.
852-
regularisation_parameter: float
853-
The main regularisation parameter to control the amount of smoothing/noise removal. Larger values lead to stronger smoothing.
854-
regularisation_iterations: int
855-
The number of iterations for regularisers (aka INNER iterations).
856-
regularisation_half_precision: bool
857-
Perform faster regularisation computation in half-precision with a very minimal sacrifice in quality.
858-
nonnegativity : bool
859-
Impose nonnegativity constraint (set to True) on the reconstructed image. Default False.
860-
gpu_id : int
861-
A GPU device index to perform operation on.
862-
863-
Returns
864-
-------
865-
cp.ndarray
866-
The ADMM reconstructed volume as a CuPy array.
867-
"""
868-
if initialisation not in ["FBP", "CGLS", None]:
869-
raise ValueError(
870-
"The acceptable values for initialisation are 'FBP','CGLS' and None"
871-
)
872-
873-
if initialisation is not None:
874-
if detector_pad == True:
875-
detector_pad = __estimate_detectorHoriz_padding(data.shape[2])
876-
877-
if detector_pad > 0:
878-
# if detector_pad is not zero we need to reconstruct the image on the recon+2*detector_pad size
879-
recon_size = data.shape[2] + 2 * detector_pad
880-
881-
if initialisation == "FBP":
780+
elif initialisation == "SIRT":
882781
initialisation_vol = cp.require(
883782
cp.swapaxes(
884-
FBP3d_tomobar(
783+
SIRT3d_tomobar(
885784
data,
886785
angles=angles,
887786
center=center,
888787
detector_pad=detector_pad,
889788
recon_size=recon_size,
890789
recon_mask_radius=recon_mask_radius,
891-
),
892-
0,
893-
1,
894-
),
895-
requirements="C",
896-
)
897-
elif initialisation == "CGLS":
898-
initialisation_vol = cp.require(
899-
cp.swapaxes(
900-
CGLS3d_tomobar(
901-
data,
902-
angles=angles,
903-
center=center,
904-
detector_pad=detector_pad,
905-
recon_size=recon_size,
906-
recon_mask_radius=recon_mask_radius,
907-
iterations=15,
790+
iterations=150,
791+
nonnegativity=True,
908792
),
909793
0,
910794
1,
@@ -915,12 +799,18 @@ def _instantiate_direct_recon_class(
915799
initialisation_vol = None
916800

917801
RecToolsCP = _instantiate_iterative_recon_class(
918-
data, angles, center, detector_pad, recon_size, gpu_id, datafidelity="LS"
802+
data,
803+
angles,
804+
center,
805+
detector_pad,
806+
recon_size,
807+
subsets_number,
808+
gpu_id,
919809
)
920810

921811
_data_ = {
922-
"projection_norm_data": data,
923-
"OS_number": subsets_number,
812+
"projection_data": data,
813+
"data_fidelity": data_fidelity,
924814
"data_axes_labels_order": input_data_axis_labels,
925815
}
926816

@@ -1044,8 +934,8 @@ def _instantiate_iterative_recon_class(
1044934
center: Optional[float] = None,
1045935
detector_pad: Union[bool, int] = False,
1046936
recon_size: Optional[int] = None,
937+
OS_number: int = 1,
1047938
gpu_id: int = 0,
1048-
datafidelity: str = "LS",
1049939
) -> Type:
1050940
"""instantiate ToMoBAR's iterative recon class
1051941
@@ -1055,8 +945,8 @@ def _instantiate_iterative_recon_class(
1055945
center (Optional[float], optional): center of recon. Defaults to None.
1056946
detector_pad : (Union[bool, int]) : Detector width padding. Defaults to False.
1057947
recon_size (Optional[int], optional): recon_size. Defaults to None.
1058-
datafidelity (str, optional): Data fidelity
1059-
gpu_id (int, optional): gpu ID. Defaults to 0.
948+
OS_number (int): The number of ordered subsets, Defaults to 1.
949+
gpu_id (int): gpu ID. Defaults to 0.
1060950
1061951
Returns:
1062952
Type[RecToolsIRCuPy]: an instance of the iterative class
@@ -1072,14 +962,14 @@ def _instantiate_iterative_recon_class(
1072962
RecToolsCP = RecToolsIRCuPy(
1073963
DetectorsDimH=data.shape[2], # Horizontal detector dimension
1074964
DetectorsDimH_pad=detector_pad, # padding for horizontal detector
1075-
DetectorsDimV=data.shape[1], # Vertical detector dimension (3D case)
965+
DetectorsDimV=data.shape[1], # Vertical detector dimension
1076966
CenterRotOffset=data.shape[2] / 2
1077967
- center
1078968
- 0.5, # Center of Rotation scalar or a vector
1079969
AnglesVec=-angles, # A vector of projection angles in radians
1080970
ObjSize=recon_size, # Reconstructed object dimensions (scalar)
1081-
datafidelity=datafidelity,
1082-
device_projector=gpu_id,
971+
device_projector=gpu_id, # device number
972+
OS_number=OS_number, # the number of ordered subsets
1083973
)
1084974
return RecToolsCP
1085975

0 commit comments

Comments
 (0)