Normalize operation (#150)

CUDA based normalization operation integrated into cucim. Normalization supports following types. 
1- simple range based normalization
2- Atangent based normalization

closes #151

Author: shekhardw

python/cucim/src/cucim/core/operations/expose/tests/test_expose.py

@@ -1,5 +1,6 @@
 from cucim.core.operations.expose.transform import (color_jitter, image_flip,
                                                     image_rotate_90,
+                                                    normalize_data,
                                                     rand_image_flip,
                                                     rand_image_rotate_90,
                                                     rand_zoom,
@@ -11,6 +12,7 @@ def test_exposed_transforms():
     assert image_flip is not None
     assert image_rotate_90 is not None
     assert scale_intensity_range is not None
+    assert normalize_data is not None
     assert zoom is not None
     assert rand_zoom is not None
     assert rand_image_flip is not None

python/cucim/src/cucim/core/operations/expose/transform.py

@@ -13,8 +13,8 @@
 # limitations under the License.
 
 from cucim.core.operations.color import color_jitter
-from cucim.core.operations.intensity import (rand_zoom, scale_intensity_range,
-                                             zoom)
+from cucim.core.operations.intensity import (normalize_data, rand_zoom,
+                                             scale_intensity_range, zoom)
 from cucim.core.operations.spatial import (image_flip, image_rotate_90,
                                            rand_image_flip,
                                            rand_image_rotate_90)

python/cucim/src/cucim/core/operations/intensity/__init__.py

@@ -1,7 +1,9 @@
+from .normalize import normalize_data
 from .scaling import scale_intensity_range
 from .zoom import rand_zoom, zoom
 
 __all__ = [
+    "normalize_data",
     "scale_intensity_range",
     "zoom", 
     "rand_zoom"

python/cucim/src/cucim/core/operations/intensity/kernel/cuda_kernel_source.py

@@ -14,14 +14,38 @@
 
 cuda_kernel_code = r'''
 extern "C" {
+__global__ void normalize_data_by_range(float *in, float *out, \
+                                        float norm_factor, \
+                                        float min_value, \
+                                        int total_size)
+{
+  const unsigned int j = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if( j < total_size ) {
+    out[j] = norm_factor * (in[j] - min_value);
+  }
+}
+
+__global__ void normalize_data_by_atan(float *in, float *out, \
+                                       float norm_factor, \
+                                       float min_value, \
+                                       int total_size)
+{
+  const unsigned int j = blockIdx.x * blockDim.x + threadIdx.x;
+
+  if( j < total_size ) {
+    out[j] = norm_factor * atan(in[j] - min_value);
+  }
+}
+
 __global__ void scaleVolume(float* image, float* output, \
                             float x, float y, float bmin, \
                             float bmax, int W)
 {
-    const unsigned int j = blockIdx.x * blockDim.x + threadIdx.x;
-    if(j < W) {
-        output[j] = fmaxf(fminf(image[j] * x - y, bmax), bmin);
-    }
+  const unsigned int j = blockIdx.x * blockDim.x + threadIdx.x;
+  if(j < W) {
+    output[j] = fmaxf(fminf(image[j] * x - y, bmax), bmin);
+  }
 }
 
 __global__ void zoom_in_kernel(float *input_tensor, float *output_tensor, \

python/cucim/src/cucim/core/operations/intensity/normalize.py

@@ -0,0 +1,129 @@
+# Copyright (c) 2021, NVIDIA CORPORATION.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# 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.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from typing import Any
+
+import cupy
+import numpy as np
+
+from .kernel.cuda_kernel_source import cuda_kernel_code
+
+CUDA_KERNELS = cupy.RawModule(code=cuda_kernel_code)
+
+
+def normalize_data(
+    img: Any,
+    norm_constant: float,
+    min_value: float,
+    max_value: float,
+    type: str = 'range'
+) -> Any:
+    """
+    Apply intensity normalization to the input array.
+    Normalize intensities to the range of [0, norm_constant].
+
+    Parameters
+    ----------
+    img : channel first, cupy.ndarray or numpy.ndarray
+        Input data of shape (C, H, W). Can also batch process input of shape
+        (N, C, H, W). Can be a numpy.ndarray or cupy.ndarray.
+    norm_constant: float
+        Normalization range of the input data. [0, norm_constant]
+    min_value : float
+        Minimum intensity value in input data.
+    max_value : float
+        Maximum intensity value in input data.
+    type : {'range', 'atan'}
+        Type of normalization.
+
+    Returns
+    -------
+    out : cupy.ndarray or numpy.ndarray
+        Output data. Same dimensions and type as input.
+
+    Raises
+    ------
+    TypeError
+        If input 'img' is not cupy.ndarray or numpy.ndarray
+    ValueError
+        If input original intensity min and max are same
+    ValueError
+        If incorrect normalization type is invoked
+
+    Examples
+    --------
+    >>> import cucim.core.operations.intensity as its
+    >>> # input is channel first 3d array
+    >>> output_array = its.normalize_data(input_arr,
+                                          10, 0 , 255)
+    """
+    try:
+        if max_value - min_value == 0.0:
+            raise ValueError("Minimum and Maximum intensity \
+                              same in input data")
+
+        if type not in ['range', 'atan']:
+            raise ValueError("Incorrect normalization type. \
+                              Supported types are: \
+                                  range based: 1,\
+                                  atangent based: 2")
+
+        to_cupy = False
+
+        if isinstance(img, np.ndarray):
+            to_cupy = True
+            cupy_img = cupy.asarray(img, dtype=cupy.float32, order='C')
+        elif not isinstance(img, cupy.ndarray):
+            raise TypeError("img must be a cupy.ndarray or numpy.ndarray")
+        else:
+            cupy_img = cupy.ascontiguousarray(img)
+
+        if cupy_img.dtype != cupy.float32:
+            if cupy.can_cast(img.dtype, cupy.float32) is False:
+                raise ValueError(
+                    "Cannot safely cast type {cupy_img.dtype.name} to \
+                    'float32'"
+                )
+            else:
+                cupy_img = cupy_img.astype(cupy.float32)
+
+        normalize = CUDA_KERNELS.get_function("normalize_data_by_range")
+
+        if type == 'atan':
+            normalize = CUDA_KERNELS.get_function("normalize_data_by_atan")
+
+        value_range = max_value - min_value
+        norm_factor = norm_constant / value_range
+
+        total_size = int(np.prod(img.shape))
+        blockx = 128
+        gridx = int((total_size - 1) / blockx + 1)
+
+        result = cupy.empty(img.shape, dtype=cupy_img.dtype)
+
+        normalize((gridx, 1, 1), (blockx, 1, 1),
+                  (cupy_img, result, np.float32(norm_factor),
+                   np.float32(min_value),
+                   np.int32(total_size)))
+
+        if img.dtype != cupy.float32:
+            result = result.astype(img.dtype)
+
+        if to_cupy is True:
+            result = cupy.asnumpy(result)
+
+    except Exception:
+        raise
+
+    return result

python/cucim/src/cucim/core/operations/intensity/tests/normalized.png

@@ -0,0 +1,77 @@
+import os
+
+import cupy
+import numpy as np
+import pytest
+import skimage.data
+from PIL import Image
+
+import cucim.core.operations.intensity as its
+
+
+def get_input_arr():
+    img = skimage.data.astronaut()
+    arr = np.asarray(img)
+    arr = np.transpose(arr)
+    return arr
+
+
+def get_norm_data():
+    dirname = os.path.dirname(__file__)
+    img1 = Image.open(os.path.join(os.path.abspath(dirname),
+                                   "normalized.png"))
+    arr_o = np.asarray(img1)
+    arr_o = np.transpose(arr_o)
+    return arr_o
+
+
+def get_norm_atan_data():
+    dirname = os.path.dirname(__file__)
+    img1 = Image.open(os.path.join(os.path.abspath(dirname),
+                                   "normalized_atan.png"))
+    arr_o = np.asarray(img1)
+    arr_o = np.transpose(arr_o)
+    return arr_o
+
+
+def test_norm_param():
+    arr = get_input_arr()
+    with pytest.raises(ValueError):
+        its.normalize_data(arr, 10.0, 255, 255)
+    with pytest.raises(ValueError):
+        its.normalize_data(arr, 10.0, 0, 255, 'invalid')
+    with pytest.raises(TypeError):
+        img = Image.fromarray(arr.T, 'RGB')
+        its.normalize_data(img, 10.0, 0, 255)
+
+
+def test_norm_numpy_input():
+    arr = get_input_arr()
+    norm_arr = get_norm_data()
+    output = its.normalize_data(arr, 10.0, 0, 255)
+    assert np.allclose(output, norm_arr)
+
+    norm_atan_arr = get_norm_atan_data()
+    output = its.normalize_data(arr, 10000, 0, 255, 'atan')
+    assert np.allclose(output, norm_atan_arr)
+
+
+def test_norm_cupy_input():
+    arr = get_input_arr()
+    norm_arr = get_norm_data()
+    cupy_arr = cupy.asarray(arr)
+    cupy_output = its.normalize_data(cupy_arr, 10.0, 0, 255)
+    np_output = cupy.asnumpy(cupy_output)
+    assert np.allclose(np_output, norm_arr)
+
+
+def test_norm_batchinput():
+    arr = get_input_arr()
+    norm_arr = get_norm_data()
+    arr_batch = np.stack((arr,) * 8, axis=0)
+    output = its.normalize_data(arr_batch, 10.0, 0, 255)
+
+    assert output.shape[0] == 8
+
+    for i in range(output.shape[0]):
+        assert np.allclose(output[i], norm_arr)