diff --git a/gratopy/fanbeam.cl b/gratopy/fanbeam.cl index e435c70..15098e7 100644 --- a/gratopy/fanbeam.cl +++ b/gratopy/fanbeam.cl @@ -566,3 +566,326 @@ __kernel void single_line_fanbeam_\my_variable_type_\order1\order2( sino[pos_sino_\order1(s, a, z, Ns, Na, Nz)] = acc * sdpd[s] / delta_xi * delta_x; } + +// Helper function for ray-driven fanbeam transforms +real fanbeam_ray_weightfkt_\my_variable_type_\order1\order2(real t,real kappa, + real s_under,real s_upper,real difference){ + + real rhs=0; + real epsilon =0.0001; + if ( fabs(difference)<(epsilon*s_under) && (fabs(t) (real)0.) { + acc += weight * img[0]; + } + // update image to next position + img += stride_x; + } + } + // assign value to sinogram + sino[pos_sino_\order1(s, a, z, Ns, Na, Nz)] = + acc * delta_x; +} + +// ray-driven Fanbeam backprojection +// Computes the backprojection projection in fanbeam geometry for a given image. +// the \my_variable_type_\order1\order2 suffix sets the kernel to the suitable +// precision, contiguity of the image, contiguity of sinogram. +// Input: +// img: Pointer to array representing image (to be computed) of +// dimensions Nx times Ny times Nz (img_shape=Nx times Ny) +// sino: Pointer to array representing sinogram (to be +// transformed) with detector-dimension times angle-dimension +// times z dimension +// ofs: Buffer containing geometric informations concerning +// the projection directions (angular information), +// The first two entries are (xd,yd) in direction along the +// detector line with length delta_xi. +// Third and fourth entries (qx,qy) from origin to source with +// length RE. +// Fifth and sixth entries (dx0,dy0) from origin to center of +// detector line (orthogonal projection of origin onto +// detector line). +// sdpd: Buffer containing the values associated with sqrt(xi^2+R^2) +// as weighting +// Geometry_information: Contains various geometric quantities +// relevant for the computation, more precisely +// 0. R/delta_x, 1. RE/delta_x, 2. delta_xi/delta_x, +// 3. x_mid (in image_pixels), 4. y_mid (in image_pixels), +// 5. xi_midpoint (in detector_pixels), 6. Nx, 7. Ny, 8. Nxi, +// 9. Nphi, 10. delta_x +// Output: +// values inside img are altered to represent the computed +// fanbeam backprojection +__kernel void fanbeam_ray_ad_\my_variable_type_\order1\order2( + __global real *img, __global real *sino, __constant real8 *ofs, + __constant real *sdpd, __constant real *Geometryinformation) { + // Extract geometric information + size_t Nx = get_global_size(0); + size_t Ny = get_global_size(1); + size_t Nz = get_global_size(2); + int Ns = Geometryinformation[8]; + int Na = Geometryinformation[9]; + + // Extract current position + size_t xx = get_global_id(0); + size_t yy = get_global_id(1); + size_t z = get_global_id(2); + + // Midpoints of geometry + real2 midpoint = (real2)(Geometryinformation[3], Geometryinformation[4]); + real midpoint_det = Geometryinformation[5]; + + // Relevant distances + real R = Geometryinformation[0]; + real delta_xi = Geometryinformation[2]; + real delta_x = Geometryinformation[10]; + + // Pixel center relative to the image midpoint + real2 P = (real2)(xx,yy) - midpoint; + + // Accumulation variable + real acc = (real)0.; + + // Shift sinogram to suitable z-position + sino += pos_sino_\order2(0, 0, z, Ns, Na, Nz); + + // precompute scaling parameters + real R_sqr_inv = (real)1. / (R * R); + real delta_xi_sqr_inv = (real)1. / (delta_xi * delta_xi); + + // for loop through all angles + for (int a = 0; a < Na; a++) { + // Geometric information associated with a.th angle + real8 o = ofs[a]; + //(xd,yd) ... vector along the detector-line with length delta_xi. + real2 dl = (real2)(o.s0, o.s1); + //(qx,qy) ... vector from origin to source (with length RE). + real2 q = (real2)(o.s2, o.s3); + //(dx0,dy0) ... vector from origin orthogonally projected onto + //detector-line.(with length R-RE) + real2 d0 = (real2)(o.s4, o.s5); + + // Delta_Phi angular resolution + real Delta_Phi = o.s6; + + real2 dd = (d0 - q) * R_sqr_inv; + real2 d = dl * delta_xi_sqr_inv; + + // define the corners + real2 Pmm = P + (real2)(-(real)0.5, -(real)0.5); + real2 Pmp = P + (real2)(-(real)0.5, (real)0.5); + real2 Ppm = P + (real2)((real)0.5, -(real)0.5); + real2 Ppp = P + (real2)((real)0.5, (real)0.5); + + //project corners onto the detector + real xi_mm = dot(d, Pmm - q) / dot(dd, Pmm - q) + midpoint_det; + real xi_mp = dot(d, Pmp - q) / dot(dd, Pmp - q) + midpoint_det; + real xi_pm = dot(d, Ppm - q) / dot(dd, Ppm - q) + midpoint_det; + real xi_pp = dot(d, Ppp - q) / dot(dd, Ppp - q) + midpoint_det; + + real xi_low = min(min(xi_mm, xi_mp), min(xi_pm, xi_pp)); + real xi_high = max(max(xi_mm, xi_mp), max(xi_pm, xi_pp)); + + int s_low = (int)floor(xi_low - (real)1.); + int s_high = (int)ceil(xi_high + (real)1.); + + s_low = max(s_low, 0); + s_high = min(s_high, Ns - 1); + + real acc_local = (real)0.; + for (int s = s_low; s <= s_high; s++) { + // Direction vector from source to the center of detector pixel s. + real2 dp = d0 + dl * ((real)s - midpoint_det) - q; + real norm = hypot(dp.x, dp.y); + dp = dp / norm; + + real2 ortho = (real2)(-dp.y, dp.x); + real abs_ortho_x = fabs(ortho.x); + real abs_ortho_y = fabs(ortho.y); + real s_under = fabs(abs_ortho_x - abs_ortho_y) * (real)0.5 * delta_x; + real s_upper = (abs_ortho_x + abs_ortho_y) * (real)0.5 * delta_x; + real difference = s_upper - s_under; + + real zz = dot(P - q, ortho) * delta_x; + real weight = fanbeam_ray_weightfkt_\my_variable_type_\order1\order2( + zz, (real)1. / max(abs_ortho_x, abs_ortho_y), s_under, s_upper, + difference); + + acc_local += weight * sino[pos_sino_\order2(s, a, 0, Ns, Na, Nz)]; + } + + // Accumulate weighted sum (Delta_Phi accounts for angular resolution) + acc += Delta_Phi * acc_local; + } + + //update img with computed value + img[pos_img_\order1(xx, yy, z, Nx, Ny, Nz)] = acc * delta_xi; +} + diff --git a/gratopy/operator/__init__.py b/gratopy/operator/__init__.py index 10a83e9..22b2e14 100644 --- a/gratopy/operator/__init__.py +++ b/gratopy/operator/__init__.py @@ -1,3 +1,9 @@ -from .projection import Radon, Fanbeam +from .projection import ( + Fanbeam, + Radon, + RayDrivenFanbeam, + RayDrivenRadon, + StripDrivenRadon, +) from .base import IDENTITY, ZERO from .opencl import OpenCLKernelSpec diff --git a/gratopy/operator/projection.py b/gratopy/operator/projection.py index 2f3885b..0533be1 100644 --- a/gratopy/operator/projection.py +++ b/gratopy/operator/projection.py @@ -40,8 +40,7 @@ from types import SimpleNamespace from typing import Any -from gratopy.gratopy import radon_struct -from gratopy.operator.base import Operator +from gratopy.gratopy import fanbeam_struct, radon_struct from gratopy.operator.opencl import OpenCLKernelSpec, _OpenCLOperator from gratopy.utilities import ( Angles, @@ -58,7 +57,6 @@ valid_image_given_detector_halfcircle, ) - class Radon(_OpenCLOperator): """Parallel-beam Radon transform operator. @@ -444,7 +442,238 @@ def apply_to( ) -class Fanbeam(Operator): +"""Fanbeam-beam Radon transform operator. + -without the substitute placeholder-""" + +class Fanbeam(_OpenCLOperator): + def __init__( + self, + source_distances: float | tuple[float, float], + image_domain: int | tuple[int, int] | ImageDomain, + angles: Angles | int, + detectors: Detectors | int | None = None, + adjoint: bool = False, + kernel_spec: OpenCLKernelSpec | None = None, + ): + if not isinstance(image_domain, ImageDomain): + image_domain = ImageDomain(size=image_domain, extent=2.0) + + if not isinstance(angles, Angles): + angles = Angles.uniform(number=angles) + + if not isinstance(detectors, Detectors): + if detectors is None: + detectors = int(np.ceil(np.hypot(*image_domain.size))) + detector_extent = image_domain.extent + if isinstance(detector_extent, ExtentPlaceholder): + detector_extent = ExtentPlaceholder.FULL + detectors = Detectors(number=detectors, extent=detector_extent) + + if isinstance(source_distances, tuple): + source_detector_distance, source_origin_distance = source_distances + else: + source_detector_distance = source_distances + source_origin_distance = source_distances / 2 + + state = { + "source_detector_distance": source_detector_distance, + "source_origin_distance": source_origin_distance, + "image_domain": image_domain, + "angles": angles, + "detectors": detectors, + "adjoint": adjoint, + } + super().__init__(name="Fanbeam", state=state, kernel_spec=kernel_spec) + + self.substitute_placeholder() + self.projection_settings: SimpleNamespace | None = None + self._host_struct: dict[str, Any] | None = None + self._device_struct: dict[tuple[cl.Context, np.dtype], dict[str, Any]] = {} + + image_shape = self.image_domain.size + sinogram_shape = (self.detectors.number, len(self.angles)) + + if self.adjoint: + self.input_shape = sinogram_shape + self.output_shape = image_shape + else: + self.input_shape = image_shape + self.output_shape = sinogram_shape + + def _default_kernel_spec(self) -> OpenCLKernelSpec: + kernel_path = Path(__file__).resolve().parent.parent / "fanbeam.cl" + return OpenCLKernelSpec.from_path(kernel_path, base_name="fanbeam") + + def __getstate__(self) -> dict[str, Any]: + """Return pickle/deepcopy state without live OpenCL runtime objects.""" + state = super().__getstate__() + state["projection_settings"] = None + state["_host_struct"] = None + state["_device_struct"] = {} + return state + + @property + def source_detector_distance(self) -> float: + return self.state["source_detector_distance"] + + @property + def source_origin_distance(self) -> float: + return self.state["source_origin_distance"] + + @property + def image_domain(self) -> ImageDomain: + return self.state["image_domain"] + + @property + def angles(self) -> Angles: + return self.state["angles"] + + @property + def detectors(self) -> Detectors: + return self.state["detectors"] + + @property + def adjoint(self) -> bool: + return self.state["adjoint"] + + @property + def T(self) -> "Fanbeam": + operator_copy = copy(self) + operator_copy.state = copy(self.state) + operator_copy.state["adjoint"] = not self.state["adjoint"] + operator_copy.input_shape, operator_copy.output_shape = ( + operator_copy.output_shape, + operator_copy.input_shape, + ) + return operator_copy + + def _repr_name_(self) -> str: + return "Fanbeam.T" if self.adjoint else "Fanbeam" + + def substitute_placeholder(self) -> None: + """(NOT IMPLEMENTED)""" + if isinstance(self.image_domain.extent, ExtentPlaceholder) or isinstance( + self.detectors.extent, ExtentPlaceholder + ): + raise NotImplementedError( + " not implemented" + ) + + def _ensure_host_struct(self, queue: cl.CommandQueue) -> None: + if self._host_struct is not None: + return + + self._host_struct = fanbeam_struct( + queue=queue, + img_shape=self.image_domain.size, + angles=self.angles.angles, + detector_width=float(self.detectors.extent), + source_detector_dist=float(self.source_detector_distance), + source_origin_dist=float(self.source_origin_distance), + angle_weights=self.angles.weights, + n_detectors=self.detectors.number, + detector_shift=self.detectors.center, + image_width=float(self.image_domain.extent), + midpoint_shift=self.image_domain.center, + reverse_detector=self.detectors.reversed, + ) + + def _ensure_device_struct( + self, + queue: cl.CommandQueue, + dtype: npt.DTypeLike, + ) -> dict[str, Any]: + self._ensure_host_struct(queue) + dtype = np.dtype(dtype) + cache_key = (queue.context, dtype) + if cache_key in self._device_struct: + return self._device_struct[cache_key] + + assert self._host_struct is not None + ofs = self._host_struct["ofs_dict"][dtype] + sdpd = self._host_struct["sdpd_dict"][dtype] + geometry = self._host_struct["geo_dict"][dtype] + angle_weights = self._host_struct["angle_diff_dict"][dtype] + + ofs_buf = cl.Buffer(queue.context, cl.mem_flags.READ_ONLY, ofs.nbytes) + sdpd_buf = cl.Buffer(queue.context, cl.mem_flags.READ_ONLY, sdpd.nbytes) + geometry_buf = cl.Buffer(queue.context, cl.mem_flags.READ_ONLY, geometry.nbytes) + angle_weights_buf = cl.Buffer( + queue.context, + cl.mem_flags.READ_ONLY, + angle_weights.nbytes, + ) + + cl.enqueue_copy(queue, ofs_buf, ofs.data).wait() + cl.enqueue_copy(queue, sdpd_buf, sdpd.data).wait() + cl.enqueue_copy(queue, geometry_buf, geometry.data).wait() + cl.enqueue_copy(queue, angle_weights_buf, angle_weights.data).wait() + + device_struct = { + "ofs": ofs_buf, + "sdpd": sdpd_buf, + "geometry": geometry_buf, + "angle_weights": angle_weights_buf, + } + self._device_struct[cache_key] = device_struct + return device_struct + + def _kernel_arguments( + self, + output: clarray.Array, + argument: clarray.Array, + queue: cl.CommandQueue, + ) -> tuple[Any, ...]: + device_struct = self._ensure_device_struct(queue, argument.dtype) + return device_struct["ofs"], device_struct["sdpd"], device_struct["geometry"] + + def apply_to( + self, + argument: npt.ArrayLike | clarray.Array, + output: clarray.Array | None = None, + queue: cl.CommandQueue | None = None, + return_event: bool = False, + ) -> clarray.Array | tuple[clarray.Array, list[cl.Event]]: + queue = self._infer_queue(argument=argument, output=output, queue=queue) + self.projection_settings = SimpleNamespace(queue=queue) + return super().apply_to( + argument, + output=output, + queue=queue, + return_event=return_event, + ) + + +"""Derived operators: ray-driven Radon and Fanbeam and strip Radon""" + +class RayDrivenRadon(Radon): + """Ray-driven Radon transform operator.""" + + def __init__( + self, + image_domain: int | tuple[int, int] | ImageDomain, + angles: Angles | int, + detectors: Detectors | int | None = None, + adjoint: bool = False, + kernel_spec: OpenCLKernelSpec | None = None, + ): + super().__init__( + image_domain=image_domain, + angles=angles, + detectors=detectors, + adjoint=adjoint, + kernel_spec=kernel_spec, + ) + self.name = "RayDrivenRadon" + + def _default_kernel_spec(self) -> OpenCLKernelSpec: + kernel_path = Path(__file__).resolve().parent.parent / "radon.cl" + return OpenCLKernelSpec.from_path(kernel_path, base_name="radon_ray") + + +class RayDrivenFanbeam(Fanbeam): + """Ray-driven fanbeam transform operator.""" + def __init__( self, source_distances: float | tuple[float, float], @@ -452,8 +681,46 @@ def __init__( angles: Angles | int, detectors: Detectors | int | None = None, adjoint: bool = False, + kernel_spec: OpenCLKernelSpec | None = None, ): - super().__init__(name="Fanbeam") + super().__init__( + source_distances=source_distances, + image_domain=image_domain, + angles=angles, + detectors=detectors, + adjoint=adjoint, + kernel_spec=kernel_spec, + ) + self.name = "RayDrivenFanbeam" + + def _default_kernel_spec(self) -> OpenCLKernelSpec: + kernel_path = Path(__file__).resolve().parent.parent / "fanbeam.cl" + return OpenCLKernelSpec.from_path(kernel_path, base_name="fanbeam_ray") + + +class StripDrivenRadon(Radon): + """Strip-driven parallel-beam Radon transform operator.""" + + def __init__( + self, + image_domain: int | tuple[int, int] | ImageDomain, + angles: Angles | int, + detectors: Detectors | int | None = None, + adjoint: bool = False, + kernel_spec: OpenCLKernelSpec | None = None, + ): + super().__init__( + image_domain=image_domain, + angles=angles, + detectors=detectors, + adjoint=adjoint, + kernel_spec=kernel_spec, + ) + self.name = "StripDrivenRadon" + + def _default_kernel_spec(self) -> OpenCLKernelSpec: + kernel_path = Path(__file__).resolve().parent.parent / "radon.cl" + return OpenCLKernelSpec.from_path(kernel_path, base_name="radon_strip") # R, R_E parameters for fanbeam tranform: pass as _one_ argument, tuple diff --git a/gratopy/radon.cl b/gratopy/radon.cl index 6f3bc25..84657c3 100644 --- a/gratopy/radon.cl +++ b/gratopy/radon.cl @@ -339,3 +339,554 @@ __kernel void single_line_radon_\my_variable_type_\order1\order2( sino[pos_sino_\order1(s, a, z, Ns, Na, Nz)] = acc * delta_x * delta_x / delta_xi; } + +// Helper function for ray-driven transforms +real ray_weightfkt_\my_variable_type_\order1\order2(real t,real kappa, + real s_under,real s_upper,real difference){ + + real rhs=0; + real epsilon =0.0001; + if ( fabs(difference)<(epsilon*s_under) && (fabs(t) (real)0.) { + acc += weight * img[0]; + } + // update image to next position + img += stride_x; + } + } + // assign value to sinogram + sino[pos_sino_\order1(s, a, z, Ns, Na, Nz)] = + acc * delta_x; +} + +// ray-driven Radon backprojection +// Computes the backprojection projection in parallel beam geometry for a given +// image. the \my_variable_type_\order1\order2 suffix sets the kernel to the +// suitable precision, contiguity of the sinogram, contiguity of image. +// Input: +// img: Pointer to array representing image (to be computed) +// of dimensions Nx times Ny times Nz +// (img_shape=Nx times Ny) +// sino: Pointer to array representing sinogram (to be +// transformed) with detector-dimension times +// angle-dimension times z dimension +// ofs: Buffer containing geometric informations concerning the +// projection-directions (angular information) +// Entries are cos, sin, offset and 1/max(|cos|,|sino|) +// Geometryinformation: Contains various geometric +// information +// [delta_x, delta_xi, Nx, Ny, Ns, Na] +// Output: +// values inside img are altered to represent the computed +// Radon backprojection +__kernel void radon_ray_ad_\my_variable_type_\order1\order2( + __global real *img, __global real *sino, __constant real8 *ofs, + __constant real *Geometryinformation) { + // Extract dimensions + size_t Nx = get_global_size(0); + size_t Ny = get_global_size(1); + size_t Nz = get_global_size(2); + const int Ns = Geometryinformation[4]; + const int Na = Geometryinformation[5]; + + // Extruct current position + size_t x = get_global_id(0); + size_t y = get_global_id(1); + size_t z = get_global_id(2); + + + // Extract scales + const float delta_x = Geometryinformation[0]; + const float delta_xi = Geometryinformation[1]; + + + // Accumulation variable + real acc = (real)0.; + real2 c = (real2)(x, y); + + // shift sinogram to correct z-dimension (as this will remain fixed), + // particularly relevant for "F" contiguity of image + sino += pos_sino_\order2(0, 0, z, Ns, Na, Nz); + + // Integrate with respect to angular dimension + for (int a = 0; a < Na; a++) { + // Extract angular dimensions + real8 o = ofs[a]; + + real s_under = fabs ( fabs(o.x) - fabs(o.y) )/2.; + real s_upper = fabs ( fabs(o.x) + fabs(o.y) )/2.; + real difference = s_upper - s_under; + + real Delta_phi = o.s4; // angle_width asociated to the angle + + // compute detector position associated to (x,y) and phi=a + real s = dot(c, o.s01) + o.s2; + + // compute adjacent detector positions + int sm = floor(s); + //int sp = sm + 1; + + int s_low = floor(s-s_upper*1.01); + int s_high = ceil(s+s_upper*1.01); + + s_low=max(s_low,0); + s_high=min(s_high,Ns-1); + + real acc_local=0; + for (int p=s_low;p<=s_high;p++) + { + real zz = s-p; + real weight = ray_weightfkt_\my_variable_type_\order1\order2(zz,fabs(o.w)*delta_x/delta_xi,s_under,s_upper,difference); + acc_local += weight * sino[pos_sino_\order2(p, a, 0, Ns, Na, Nz)]; + } + + // accumulate weigthed sum (Delta_Phi weight due to angular resolution) + acc += Delta_phi * acc_local; + } + + // Assign value to img + img[pos_img_\order1(x, y, z, Nx, Ny, Nz)] = acc*delta_xi/delta_x; +} + + +// Strip Radon Helper functions +real proj_on_interval_\my_variable_type_\order1\order2(real x,real y,real z) +{//Projects the value of z onto the interval [x,y] +return max(x,min(y,z)); +} + + +real strip_central_primative_\my_variable_type_\order1\order2(real t,real kappa, real delta_x) +{//calculates the primative of the strip weight for the central interval + + +return 1/delta_x * kappa * t ; +} + +real strip_left_primative_\my_variable_type_\order1\order2(real t,real kappa, real s_upper, real delta_x,real difference,real inv_difference) +{//calculates the primative of the strip weight for the left interval + +if (fabs(difference)<0.0001 * s_upper) +{ +return 0; +} +return 1/delta_x * (t*s_upper + t*t/2)*inv_difference*kappa ; +} + +real strip_right_primative_\my_variable_type_\order1\order2(real t,real kappa, real s_upper, real delta_x,real difference,real inv_difference) +{//calculates the primative of the strip weight for the right interval + + +if (fabs(difference)<0.0001 * s_upper) +{ +return 0; +} + +return 1/delta_x * (t*s_upper - t*t/2)*inv_difference*kappa ; +} + + +real strip_weightfkt_\my_variable_type_\order1\order2(real t,real kappa, real s_under, real s_upper, real difference,real inv_difference,real delta_x) +{//calculates weight function for the strip + +real t_central = proj_on_interval_\my_variable_type_\order1\order2(-s_under,s_under, t); +real t_left = proj_on_interval_\my_variable_type_\order1\order2(-s_upper,-s_under, t); +real t_right = proj_on_interval_\my_variable_type_\order1\order2(s_under,s_upper, t); + +real val_central = strip_central_primative_\my_variable_type_\order1\order2(t_central, kappa, delta_x); +real val_left = strip_left_primative_\my_variable_type_\order1\order2(t_left,kappa, s_upper, delta_x, difference, inv_difference); +real val_right = strip_right_primative_\my_variable_type_\order1\order2(t_right,kappa, s_upper, delta_x, difference, inv_difference); + +return val_central+val_left+val_right; +} + + +// Strip-driven Radon Transform +// Computes the forward projection in parallel beam geometry for a given image. +// the \my_variable_type_\order1\order2 suffix sets the kernel to the suitable +// precision, contiguity of the sinogram, contiguity of image. +// Input: +// sino: Pointer to array representing sinogram (to be +// computed) with detector-dimension times +// angle-dimension times z dimension +// img: Pointer to array representing image to be transformed +// of dimensions Nx times Ny times Nz +// (img_shape=Nx times Ny) +// ofs: Buffer containing geometric informations concerning the +// projection-directions (angular information) +// Entries are cos, sin, offset and 1/max(|cos|,|sino|) +// Geometryinformation: Contains various geometric +// information +// [delta_x, delta_xi, Nx, Ny, Ns, Na] +// Output: +// values inside sino are altered to represent the computed +// Radon transform +__kernel void radon_strip_\my_variable_type_\order1\order2( + __global real *sino, __global real *img, __constant real8 *ofs, + __constant real *Geometryinformation) { + // Extract dimensions + size_t Ns = get_global_size(0); + size_t Na = get_global_size(1); + size_t Nz = get_global_size(2); + const int Nx = Geometryinformation[2]; + const int Ny = Geometryinformation[3]; + + // Extract current position + size_t s = get_global_id(0); + size_t a = get_global_id(1); + size_t z = get_global_id(2); + + + + + // Extract scales + const float delta_x = Geometryinformation[0]; + const float delta_xi = Geometryinformation[1]; + + // hack (since otherwise s is unsigned which leads to overflow problems) + int ss = s; + + // Extract angular information + // o = (cos,sin,offset,1/max(|cos|,|sin|)) + real4 o = ofs[a].s0123; + real reverse_mask = ofs[a].s5; + real s_under = fabs ( fabs(o.x) - fabs(o.y) )/2.; + real s_upper = fabs ( fabs(o.x) + fabs(o.y) )/2.; + real difference = s_upper - s_under; + real inv_difference = 1.0/difference; + + + // Dummy variable for switching from horizontal to vertical lines + int Nxx = Nx; + int Nyy = Ny; + int horizontal = 1; + + // When line is horizontal rather than vertical, switch x and y dimensions + if (reverse_mask != (real)0.) { + horizontal = 0; + o.xy = (real2)(o.y, o.x); + + Nxx = Ny; + Nyy = Nx; + } + + // accumulation variable + real acc = (real)0.; + + // shift image to correct z-dimension (as this will remain fixed), + // particularly relevant for "F" contiguity of image + __global real *img0 = img + pos_img_\order2(0, 0, z, Nx, Ny, Nz); + + // stride representing one index_step in x dimension (dependent on + // horizontal/vertical) + size_t stride_x = horizontal == 1 ? pos_img_\order2(1, 0, 0, Nx, Ny, Nz) + : pos_img_\order2(0, 1, 0, Nx, Ny, Nz); + + // for through the entire y dimension + for (int y = 0; y < Nyy; y++) { + int x_low, x_high; + + // project (0,y) onto detector + real d = y * o.y + o.z - ss; + + // compute bounds: minimal enclosing integer range (ceil/floor instead of an + // (int) cast to the enclosing integer). The +-0.5 is the strip box half-width; + // the *1.01 cushion is kept for boundary robustness. min/max replaces the swap. + real x_a = (-s_upper*1.01 - 0.5 - d) * o.w; + real x_b = ( s_upper*1.01 + 0.5 - d) * o.w; + x_low = (int)ceil (min(x_a, x_b)); + x_high = (int)floor(max(x_a, x_b)); + + // make sure x inside image dimensions + x_low = max(x_low, 0); + x_high = min(x_high, Nxx - 1); + + // shift position of image depending on horizontal/vertical + if (horizontal == 1) + img = img0 + pos_img_\order2(x_low, y, 0, Nx, Ny, Nz); + if (horizontal == 0) + img = img0 + pos_img_\order2(y, x_low, 0, Nx, Ny, Nz); + + // integration in x dimension for fixed y + for (int x = x_low; x <= x_high; x++) { + // anterpolation weight via normal distance + real zz = x * o.x + d; + + real weight=0; + real weight1 = strip_weightfkt_\my_variable_type_\order1\order2(zz-0.5,fabs(o.w)*delta_x/delta_xi,s_under,s_upper,difference,inv_difference,delta_x/delta_xi); + real weight2 = strip_weightfkt_\my_variable_type_\order1\order2(zz+0.5,fabs(o.w)*delta_x/delta_xi,s_under,s_upper,difference,inv_difference,delta_x/delta_xi); + + weight = weight2-weight1; + + if ( (a==900) && (ss==2000)) + { + // printf ("### Weight=%f,\n", weight); + } + + + //if (weight > (real)0.) { + acc += weight * img[0]; + //} + // update image to next position + img += stride_x; + } + } + // assign value to sinogram + // delta_x*delta_x/delta_xi (as in the standard/pixel kernel): the strip + // weight density collapses to |o.w| (the 1/(delta_x/delta_xi) in the strip + // primitive cancels the delta_x/delta_xi inside kappa), so the extra + // delta_x/delta_xi factor that the ray weight carries must be supplied here. + // Without it the strip projection scales as Nx/Ns and only matches ray/pixel + // on the diagonal Nx==Ns. + sino[pos_sino_\order1(s, a, z, Ns, Na, Nz)] = + acc * delta_x * delta_x / delta_xi; +} + +// Strip-driven Radon backprojection ##### (NOT IMPLEMENTED YET) +// Computes the backprojection projection in parallel beam geometry for a given +// image. the \my_variable_type_\order1\order2 suffix sets the kernel to the +// suitable precision, contiguity of the sinogram, contiguity of image. +// Input: +// img: Pointer to array representing image (to be computed) +// of dimensions Nx times Ny times Nz +// (img_shape=Nx times Ny) +// sino: Pointer to array representing sinogram (to be +// transformed) with detector-dimension times +// angle-dimension times z dimension +// ofs: Buffer containing geometric informations concerning the +// projection-directions (angular information) +// Entries are cos, sin, offset and 1/max(|cos|,|sino|) +// Geometryinformation: Contains various geometric +// information +// [delta_x, delta_xi, Nx, Ny, Ns, Na] +// Output: +// values inside img are altered to represent the computed +// Radon backprojection +__kernel void radon_strip_ad_\my_variable_type_\order1\order2( + __global real *img, __global real *sino, __constant real8 *ofs, + __constant real *Geometryinformation) { + + + + // Extract dimensions + size_t Nx = get_global_size(0); + size_t Ny = get_global_size(1); + size_t Nz = get_global_size(2); + const int Ns = Geometryinformation[4]; + const int Na = Geometryinformation[5]; + + // Extruct current position + size_t x = get_global_id(0); + size_t y = get_global_id(1); + size_t z = get_global_id(2); + + + // Extract scales + const float delta_x = Geometryinformation[0]; + const float delta_xi = Geometryinformation[1]; + + + // Accumulation variable + real acc = (real)0.; + real2 c = (real2)(x, y); + + // shift sinogram to correct z-dimension (as this will remain fixed), + // particularly relevant for "F" contiguity of image + sino += pos_sino_\order2(0, 0, z, Ns, Na, Nz); + + // Integrate with respect to angular dimension + for (int a = 0; a < Na; a++) { + // Extract angular dimensions + real8 o = ofs[a]; + + real s_under = fabs ( fabs(o.x) - fabs(o.y) )/2.; + real s_upper = fabs ( fabs(o.x) + fabs(o.y) )/2.; + real difference = s_upper - s_under; + real inv_difference = 1.0/difference; + + real Delta_phi = o.s4; // angle_width asociated to the angle + + // compute detector position associated to (x,y) and phi=a + real s = dot(c, o.s01) + o.s2; + + // compute adjacent detector positions + //int sm = floor(s); + //int sp = sm + 1; + + // minimal enclosing integer range (ceil/floor instead of floor/ceil) -> no + // extra zero-weight detector pixel per side. The +-0.5 is the strip's box + // half-width; the *1.01 cushion is kept for boundary robustness. + int s_low = ceil(s-s_upper*1.01-0.5); + int s_high = floor(s+s_upper*1.01+0.5); + + s_low=max(s_low,0); + s_high=min(s_high,Ns-1); + // if ((x==500)&&(y==500)) + // {printf ("ADSF: %d , %d, %f, %f \n", s_low,s_high,s,s_upper);} + + real acc_local=0; + for (int p=s_low;p<=s_high;p++) + { + real zz = s-p; + real weight=0; + real weight1 = strip_weightfkt_\my_variable_type_\order1\order2(zz-0.5,fabs(o.w)*delta_x/delta_xi,s_under,s_upper,difference,inv_difference,delta_x/delta_xi); + real weight2 = strip_weightfkt_\my_variable_type_\order1\order2(zz+0.5,fabs(o.w)*delta_x/delta_xi,s_under,s_upper,difference,inv_difference,delta_x/delta_xi); + + weight = weight2-weight1; + + acc_local += weight * sino[pos_sino_\order2(p, a, 0, Ns, Na, Nz)]; + } + + + // accumulate weigthed sum (Delta_Phi weight due to angular resolution) + acc += Delta_phi * acc_local; + } + + // Assign value to img + // analogous to the forward strip fix: the strip weight density collapses to + // |o.w| (missing the delta_x/delta_xi the ray weight carries), so the extra + // delta_x/delta_xi factor must be supplied here. Combined with the ray-adjoint + // scaling delta_xi/delta_x this collapses to a plain acc. Without it the strip + // backprojection of a constant scales as Nx/Ns and only matches R*(1)=pi on + // the diagonal Nx==Ns. + img[pos_img_\order1(x, y, z, Nx, Ny, Nz)] = acc; +} + diff --git a/tests/test_fanbeam.py b/tests/test_fanbeam.py index c0b03da..18aee9b 100644 --- a/tests/test_fanbeam.py +++ b/tests/test_fanbeam.py @@ -12,6 +12,7 @@ from pathlib import Path from .helpers import evaluate_control_numbers, create_phantoms +from gratopy.utilities import Detectors, ImageDomain # Plots are deactivated by default, can be activated @@ -1578,6 +1579,218 @@ def test_create_sparse_matrix(dtype): name="backprojected image", ) +# New operator basic tests written generally since fanbeam strip driven might be added in the future + +def _test_operator_projection(operator_class, name): + """ + Basic projection test for operator API fanbeam geometry variants. + + Computes forward and backprojection for two test images and repeats both + calls to estimate execution time, mirroring test_projection above. + """ + + print(name + " projection test") + + ctx = cl.create_some_context(interactive=INTERACTIVE) + queue = cl.CommandQueue(ctx) + + dtype = np.dtype("float32") + N = 1200 + img_gpu = create_phantoms(queue, N, dtype=dtype) + + angles = 360 + number_detectors = 600 + detector_width = 400.0 + source_detector_distance = 752.0 + source_origin_distance = 200.0 + + PS = gratopy.ProjectionSettings( + queue, + gratopy.FANBEAM, + img_shape=img_gpu.shape, + angles=angles, + detector_width=detector_width, + RE=source_origin_distance, + R=source_detector_distance, + n_detectors=number_detectors, + ) + projection = operator_class( + source_distances=(source_detector_distance, source_origin_distance), + image_domain=ImageDomain(size=(N, N), extent=PS.image_width), + angles=angles, + detectors=Detectors(number=number_detectors, extent=detector_width), + ) + backprojection = projection.T + + sino_gpu = clarray.zeros( + queue, projection.output_shape + (2,), dtype=dtype, order="F" + ) + backprojected_gpu = clarray.zeros( + queue, projection.input_shape + (2,), dtype=dtype, order="F" + ) + + iterations = 10 + a = time.perf_counter() + for i in range(iterations): + projection.apply_to(img_gpu, output=sino_gpu) + sino_gpu.get() + print( + "Average time required for " + + name + + " forward projection " + + f"{(time.perf_counter() - a) / iterations:.3f}" + ) + + a = time.perf_counter() + for i in range(iterations): + backprojection.apply_to(sino_gpu, output=backprojected_gpu) + backprojected_gpu.get() + print( + "Average time required for " + + name + + " backprojection " + + f"{(time.perf_counter() - a) / iterations:.3f}" + ) + + img = img_gpu.get() + sino = sino_gpu.get() + backprojected = backprojected_gpu.get() + + if PLOT: + plt.figure() + plt.imshow(np.hstack([img[:, :, 0], img[:, :, 1]]), cmap=plt.cm.gray) + plt.title(name + " original image") + plt.figure() + plt.imshow(np.hstack([sino[:, :, 0], sino[:, :, 1]]), cmap=plt.cm.gray) + plt.title(name + " sinogram") + plt.figure() + plt.imshow( + np.hstack([backprojected[:, :, 0], backprojected[:, :, 1]]), + cmap=plt.cm.gray, + ) + plt.title(name + " backprojected image") + plt.show() + + evaluate_control_numbers( + img, + (N, N, number_detectors, angles, 2), + expected_result=2949.3738, + classified="img", + name=name + " original image", + ) + + assert np.all(np.isfinite(sino)) + assert np.linalg.norm(sino) > 0 + assert np.all(np.isfinite(backprojected)) + assert np.linalg.norm(backprojected) > 0 + + +def test_projection_ray(): + _test_operator_projection( + gratopy.operator.RayDrivenFanbeam, "ray-driven Fanbeam" + ) + + +def _test_operator_adjointness(operator_class, name): + """Adjointness test for operator API fanbeam derived fanbeam geometry variants.""" + + print(name + " adjointness test") + + ctx = cl.create_some_context(interactive=INTERACTIVE) + queue = cl.CommandQueue(ctx) + + dtype = np.dtype("float32") + order = "F" + + Nx = 400 + number_detectors = 230 + angles = 360 + img_shape = (Nx, Nx) + detector_width = 83.0 + detector_shift = 0.0 + midpoint_shift = [0.0, 0.0] + source_detector_distance = 900.0 + source_origin_distance = 300.0 + + PS = gratopy.ProjectionSettings( + queue, + gratopy.FANBEAM, + img_shape, + angles, + n_detectors=number_detectors, + detector_width=detector_width, + detector_shift=detector_shift, + midpoint_shift=midpoint_shift, + R=source_detector_distance, + RE=source_origin_distance, + image_width=None, + ) + projection = operator_class( + source_distances=(source_detector_distance, source_origin_distance), + image_domain=ImageDomain(size=img_shape, extent=PS.image_width), + angles=angles, + detectors=Detectors( + number=number_detectors, + extent=detector_width, + center=detector_shift, + ), + ) + backprojection = projection.T + + sino2_gpu = clarray.zeros(queue, projection.output_shape, dtype=dtype, order=order) + img2_gpu = clarray.zeros(queue, projection.input_shape, dtype=dtype, order=order) + + Error = [] + count = 0 + eps = 0.00001 + + for i in range(100): + img1_gpu = clarray.to_device( + queue, np.require(np.random.random(projection.input_shape), dtype, order) + ) + sino1_gpu = clarray.to_device( + queue, np.require(np.random.random(projection.output_shape), dtype, order) + ) + + projection.apply_to(img1_gpu, output=sino2_gpu) + backprojection.apply_to(sino1_gpu, output=img2_gpu) + + pairing_img = clarray.vdot(img1_gpu, img2_gpu).get() * PS.delta_x**2 + pairing_sino = ( + clarray.vdot(gratopy.weight_sinogram(sino1_gpu, PS), sino2_gpu).get() + * PS.delta_s + ) + + relative_error = abs(pairing_img - pairing_sino) / min( + abs(pairing_img), abs(pairing_sino) + ) + if relative_error > eps: + count += 1 + Error.append((pairing_img, pairing_sino)) + + print( + name + + " adjointness: Number of Errors: " + + str(count) + + " out of 100 tests adjointness-errors were bigger than " + + str(eps) + ) + assert len(Error) < 10, ( + "A large number of " + + name + + " experiments for adjointness turned out negative, number of errors: " + + str(count) + + " out of 100 tests adjointness-errors were bigger than " + + str(eps) + ) + + +def test_adjointness_ray(): + _test_operator_adjointness( + gratopy.operator.RayDrivenFanbeam, "ray-driven Fanbeam" + ) + + # test if __name__ == "__main__": diff --git a/tests/test_radon.py b/tests/test_radon.py index c7a7ed6..9b45a36 100644 --- a/tests/test_radon.py +++ b/tests/test_radon.py @@ -9,6 +9,7 @@ import gratopy from .helpers import evaluate_control_numbers, create_phantoms +from gratopy.utilities import Detectors, ImageDomain # Plots are deactivated by default, can be activated # by setting 'export GRATOPY_TEST_PLOT=true' in the terminal @@ -1116,6 +1117,188 @@ def test_angle_input_variants(): ) +def _test_operator_projection(operator_class, name): + """ + Basic projection test for operator API parallel geometry Radon variants. + + Computes forward and backprojection for two test images and repeats both + calls to estimate execution time, mirroring test_projection above. + """ + + print(name + " projection test") + + ctx = cl.create_some_context(interactive=INTERACTIVE) + queue = cl.CommandQueue(ctx) + + dtype = np.dtype("float32") + N = 1200 + img_gpu = create_phantoms(queue, N, dtype=dtype) + + angles = 360 + detector_width = 4.0 + image_width = 4.0 + Ns = int(0.5 * N) + + projection = operator_class( + image_domain=ImageDomain(size=(N, N), extent=image_width), + angles=angles, + detectors=Detectors(number=Ns, extent=detector_width), + ) + backprojection = projection.T + + sino_gpu = clarray.zeros( + queue, projection.output_shape + (2,), dtype=dtype, order="F" + ) + backprojected_gpu = clarray.zeros( + queue, projection.input_shape + (2,), dtype=dtype, order="F" + ) + + M = 10 + a = time.perf_counter() + for i in range(M): + projection.apply_to(img_gpu, output=sino_gpu) + sino_gpu.get() + print( + "Average time required for " + + name + + " forward projection " + + str((time.perf_counter() - a) / M) + ) + + a = time.perf_counter() + for i in range(M): + backprojection.apply_to(sino_gpu, output=backprojected_gpu) + backprojected_gpu.get() + print( + "Average time required for " + + name + + " backprojection " + + str((time.perf_counter() - a) / M) + ) + + img = img_gpu.get() + sino = sino_gpu.get() + backprojected = backprojected_gpu.get() + + if PLOT: + plt.figure() + plt.imshow(np.hstack([img[:, :, 0], img[:, :, 1]]), cmap=plt.cm.gray) + plt.title(name + " original image") + plt.figure() + plt.imshow(np.hstack([sino[:, :, 0], sino[:, :, 1]]), cmap=plt.cm.gray) + plt.title(name + " sinogram") + plt.figure() + plt.imshow( + np.hstack([backprojected[:, :, 0], backprojected[:, :, 1]]), + cmap=plt.cm.gray, + ) + plt.title(name + " backprojected image") + plt.show() + + evaluate_control_numbers( + img, + (N, N, Ns, angles, 2), + expected_result=2949.3738, + classified="img", + name=name + " original image", + ) + + assert np.all(np.isfinite(sino)) + assert np.linalg.norm(sino) > 0 + assert np.all(np.isfinite(backprojected)) + assert np.linalg.norm(backprojected) > 0 + + +def test_projection_ray(): + _test_operator_projection(gratopy.operator.RayDrivenRadon, "ray-driven Radon") + + +def test_projection_strip(): + _test_operator_projection(gratopy.operator.StripDrivenRadon, "strip-driven Radon") + + +def _test_operator_adjointness(operator_class, name): + """Adjointness test for operator API Radon variants.""" + + print(name + " adjointness test") + + ctx = cl.create_some_context(interactive=INTERACTIVE) + queue = cl.CommandQueue(ctx) + + dtype = np.dtype("float32") + order = "F" + + Nx = 400 + number_detectors = 230 + angles = 180 + img_shape = (Nx, Nx) + + PS = gratopy.ProjectionSettings( + queue, gratopy.PARALLEL, img_shape, angles, n_detectors=number_detectors + ) + projection = operator_class( + image_domain=ImageDomain(size=img_shape, extent=PS.image_width), + angles=angles, + detectors=Detectors(number=number_detectors, extent=PS.detector_width), + ) + backprojection = projection.T + + sino2_gpu = clarray.zeros(queue, projection.output_shape, dtype=dtype, order=order) + img2_gpu = clarray.zeros(queue, projection.input_shape, dtype=dtype, order=order) + + Error = [] + count = 0 + eps = 0.00001 + + for i in range(100): + img1_gpu = clarray.to_device( + queue, np.require(np.random.random(projection.input_shape), dtype, order) + ) + sino1_gpu = clarray.to_device( + queue, np.require(np.random.random(projection.output_shape), dtype, order) + ) + + projection.apply_to(img1_gpu, output=sino2_gpu) + backprojection.apply_to(sino1_gpu, output=img2_gpu) + + pairing_img = clarray.vdot(img1_gpu, img2_gpu).get() * PS.delta_x**2 + pairing_sino = ( + clarray.vdot(gratopy.weight_sinogram(sino1_gpu, PS), sino2_gpu).get() + * PS.delta_s + ) + + relative_error = abs(pairing_img - pairing_sino) / min( + abs(pairing_img), abs(pairing_sino) + ) + if relative_error > eps: + count += 1 + Error.append((pairing_img, pairing_sino)) + + print( + name + + " adjointness: Number of Errors: " + + str(count) + + " out of 100 tests adjointness-errors were bigger than " + + str(eps) + ) + assert len(Error) < 10, ( + "A large number of " + + name + + " experiments for adjointness turned out negative, number of errors: " + + str(count) + + " out of 100 tests adjointness-errors were bigger than " + + str(eps) + ) + + +def test_adjointness_ray(): + _test_operator_adjointness(gratopy.operator.RayDrivenRadon, "ray-driven Radon") + + +def test_adjointness_strip(): + _test_operator_adjointness(gratopy.operator.StripDrivenRadon, "strip-driven Radon") + + # test if __name__ == "__main__": pass