Skip to content

Commit

Permalink
corr11
Browse files Browse the repository at this point in the history
  • Loading branch information
@dkazanc
dkazanc committed May 1, 2024
1 parent aad6294 commit 7e88035
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 58 deletions.
86 changes: 46 additions & 40 deletions httomolibgpu/misc/morph.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,23 +21,8 @@
"""Module for data type morphing functions"""

import numpy as np
cupy_run = False
try:
import cupy as xp

try:
xp.cuda.Device(0).compute_capability
cupy_run = True
except xp.cuda.runtime.CUDARuntimeError:
print("CuPy library is a major dependency for HTTomolibgpu, please install")
import numpy as xp
except ImportError:
import numpy as xp

if cupy_run:
from cupyx.scipy.interpolate import interpn
else:
from scipy.interpolate import interpn
from httomolibgpu import cupywrapper
cp = cupywrapper.cp

import nvtx
from typing import Literal
Expand All @@ -47,11 +32,10 @@
"data_resampler",
]


@nvtx.annotate()
def sino_360_to_180(
data: xp.ndarray, overlap: int = 0, rotation: Literal["left", "right"] = "left"
) -> xp.ndarray:
data: cp.ndarray, overlap: int = 0, rotation: Literal["left", "right"] = "left"
) -> cp.ndarray:
"""
Converts 0-360 degrees sinogram to a 0-180 sinogram.
If the number of projections in the input data is odd, the last projection
Expand All @@ -71,6 +55,17 @@ def sino_360_to_180(
cp.ndarray
Output 3D data.
"""
if cupywrapper.cupy_run:
return __sino_360_to_180(data, overlap, rotation)
else:
print("sino_360_to_180 won't be executed because CuPy is not installed")
return data


def __sino_360_to_180(
data: cp.ndarray, overlap: int = 0, rotation: Literal["left", "right"] = "left"
) -> cp.ndarray:

if data.ndim != 3:
raise ValueError("only 3D data is supported")

Expand All @@ -84,18 +79,18 @@ def sino_360_to_180(

n = dx // 2

out = xp.empty((n, dy, 2 * dz - overlap), dtype=data.dtype)
out = cp.empty((n, dy, 2 * dz - overlap), dtype=data.dtype)

if rotation == "left":
weights = xp.linspace(0, 1.0, overlap)
weights = cp.linspace(0, 1.0, overlap)
out[:, :, -dz + overlap :] = data[:n, :, overlap:]
out[:, :, : dz - overlap] = data[n : 2 * n, :, overlap:][:, :, ::-1]
out[:, :, dz - overlap : dz] = (
weights * data[:n, :, :overlap]
+ (weights * data[n : 2 * n, :, :overlap])[:, :, ::-1]
)
elif rotation == "right":
weights = xp.linspace(1.0, 0, overlap)
weights = cp.linspace(1.0, 0, overlap)
out[:, :, : dz - overlap] = data[:n, :, :-overlap]
out[:, :, -dz + overlap :] = data[n : 2 * n, :, :-overlap][:, :, ::-1]
out[:, :, dz - overlap : dz] = (
Expand All @@ -110,8 +105,8 @@ def sino_360_to_180(

@nvtx.annotate()
def data_resampler(
data: xp.ndarray, newshape: list, axis: int = 1, interpolation: str = "linear"
) -> xp.ndarray:
data: cp.ndarray, newshape: list, axis: int = 1, interpolation: str = "linear"
) -> cp.ndarray:
"""Down/Up-resampler of the input data implemented through interpn function.
Please note that the method will leave the specified axis
dimension unchanged, e.g. (128,128,128) -> (128,256,256) for axis = 0 and
Expand All @@ -128,31 +123,42 @@ def data_resampler(
Returns:
cp.ndarray: Up/Down-scaled 3D cupy array
"""
"""
if cupywrapper.cupy_run:
return __data_resampler(data, newshape, axis, interpolation)
else:
print("data_resampler won't be executed because CuPy is not installed")
return data

def __data_resampler(
data: cp.ndarray, newshape: list, axis: int = 1, interpolation: str = "linear"
) -> cp.ndarray:

from cupyx.scipy.interpolate import interpn

if data.ndim != 3:
raise ValueError("only 3D data is supported")

N, M, Z = xp.shape(data)
N, M, Z = cp.shape(data)

if axis == 0:
xaxis = xp.arange(M) - M / 2
yaxis = xp.arange(Z) - Z / 2
xaxis = cp.arange(M) - M / 2
yaxis = cp.arange(Z) - Z / 2
step_x = M / newshape[0]
step_y = Z / newshape[1]
scaled_data = xp.empty((N, newshape[0], newshape[1]), dtype=xp.float32)
scaled_data = cp.empty((N, newshape[0], newshape[1]), dtype=cp.float32)
elif axis == 1:
xaxis = xp.arange(N) - N / 2
yaxis = xp.arange(Z) - Z / 2
xaxis = cp.arange(N) - N / 2
yaxis = cp.arange(Z) - Z / 2
step_x = N / newshape[0]
step_y = Z / newshape[1]
scaled_data = xp.empty((newshape[0], M, newshape[1]), dtype=xp.float32)
scaled_data = cp.empty((newshape[0], M, newshape[1]), dtype=cp.float32)
elif axis == 2:
xaxis = xp.arange(N) - N / 2
yaxis = xp.arange(M) - M / 2
xaxis = cp.arange(N) - N / 2
yaxis = cp.arange(M) - M / 2
step_x = N / newshape[0]
step_y = M / newshape[1]
scaled_data = xp.empty((newshape[0], newshape[1], Z), dtype=xp.float32)
scaled_data = cp.empty((newshape[0], newshape[1], Z), dtype=cp.float32)
else:
raise ValueError("Only 0,1,2 values for axes are supported")

Expand Down Expand Up @@ -181,7 +187,7 @@ def data_resampler(
xi_size = xi.size
xi = np.rollaxis(xi, 0, 3)
xi = np.reshape(xi, [xi_size // 2, 2])
xi = xp.asarray(xi, dtype=xp.float32, order="C")
xi = cp.asarray(xi, dtype=cp.float32, order="C")

if axis == 0:
for j in range(N):
Expand All @@ -193,7 +199,7 @@ def data_resampler(
bounds_error=False,
fill_value=0.0,
)
scaled_data[j, :, :] = xp.reshape(
scaled_data[j, :, :] = cp.reshape(
res, [newshape[0], newshape[1]], order="C"
)
elif axis == 1:
Expand All @@ -207,7 +213,7 @@ def data_resampler(
bounds_error=False,
fill_value=0.0,
)
scaled_data[:, j, :] = xp.reshape(
scaled_data[:, j, :] = cp.reshape(
res, [newshape[0], newshape[1]], order="C"
)
else:
Expand All @@ -220,7 +226,7 @@ def data_resampler(
bounds_error=False,
fill_value=0.0,
)
scaled_data[:, :, j] = xp.reshape(
scaled_data[:, :, j] = cp.reshape(
res, [newshape[0], newshape[1]], order="C"
)

Expand Down
70 changes: 52 additions & 18 deletions httomolibgpu/misc/rescale.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either ecpress or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ---------------------------------------------------------------------------
Expand All @@ -20,28 +20,62 @@
# ---------------------------------------------------------------------------

import numpy as np
try:
import cupy as xp

try:
xp.cuda.Device(0).compute_capability
except xp.cuda.runtime.CUDARuntimeError:
print("CuPy library is a major dependency for HTTomolibgpu, please install")
import numpy as xp
except ImportError:
import numpy as xp
from httomolibgpu import cupywrapper
cp = cupywrapper.cp

import nvtx
from typing import Literal, Optional, Tuple, Union

__all__ = [
"rescale_to_int",
]


@nvtx.annotate()
def rescale_to_int(
data: xp.ndarray,
data: cp.ndarray,
perc_range_min: float = 0.0,
perc_range_max: float = 100.0,
bits: Literal[8, 16, 32] = 8,
glob_stats: Optional[Tuple[float, float, float, int]] = None,
):
"""
Rescales the data and converts it fit into the range of an unsigned integer type
with the given number of bits.
Parameters
----------
data : cp.ndarray
Required input data array, on GPU
perc_range_min: float, optional
The lower cutoff point in the input data, in percent of the data range (defaults to 0).
The lower bound is computed as min + perc_range_min/100*(max-min)
perc_range_max: float, optional
The upper cutoff point in the input data, in percent of the data range (defaults to 100).
The upper bound is computed as min + perc_range_max/100*(max-min)
bits: Literal[8, 16, 32], optional
The number of bits in the output integer range (defaults to 8).
Allowed values are:
- 8 -> uint8
- 16 -> uint16
- 32 -> uint32
glob_stats: tuple, optional
Global statistics of the full dataset (beyond the data passed into this call).
It's a tuple with (min, max, sum, num_items). If not given, the min/max is
computed from the given data.
Returns
-------
cp.ndarray
The original data, clipped to the range specified with the perc_range_min and
perc_range_max, and scaled to the full range of the output integer type
"""
if cupywrapper.cupy_run:
return __rescale_to_int(data, perc_range_min, perc_range_max, bits, glob_stats)
else:
print("rescale_to_int won't be executed because CuPy is not installed")
return data

def __rescale_to_int(
data: cp.ndarray,
perc_range_min: float = 0.0,
perc_range_max: float = 100.0,
bits: Literal[8, 16, 32] = 8,
Expand Down Expand Up @@ -91,8 +125,8 @@ def rescale_to_int(
output_max = np.iinfo(output_dtype).max

if not isinstance(glob_stats, tuple):
min_value = float(xp.min(data))
max_value = float(xp.max(data))
min_value = float(cp.min(data))
max_value = float(cp.max(data))
else:
min_value = glob_stats[0]
max_value = glob_stats[1]
Expand All @@ -103,8 +137,8 @@ def rescale_to_int(

factor = (output_max - output_min) / (input_max - input_min)

res = xp.empty(data.shape, dtype=output_dtype)
rescale_kernel = xp.ElementwiseKernel(
res = cp.empty(data.shape, dtype=output_dtype)
rescale_kernel = cp.ElementwiseKernel(
"T x, raw T input_min, raw T input_max, raw T factor",
"O out",
"""
Expand Down

0 comments on commit 7e88035

Please sign in to comment.