Repository: Eve-ning/glcm-cupy (language: python)

Issue:
Combine GLCM class to glcm and GLCM.run variants



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index d172fbe..4b3c62d 100644
--- a/README.md
+++ b/README.md
@@ -4,12 +4,147 @@ This is a CuPy reimplementation of my `glcmbin5` [**on my other repository**](ht
 
 This directly utilizes CUDA to speed up the processing of GLCM.
 
-## Notes
+# Installation
+
+First, you need to install this
+
+```shell
+pip install glcm-cupy
+```
+
+Then, you need **CuPy**.
+You need to [install CuPy manually](https://docs.cupy.dev/en/stable/install.html), 
+as it's dependent on the version of CUDA you have.
+
+I recommend using `conda-forge` as it worked for me :)
+
+For CUDA `11.6`, we use
+```shell
+conda install -c conda-forge cupy cudatoolkit=11.6
+```
+
+Replace the version you have on the arg.
+
+```shell
+conda install -c conda-forge cupy cudatoolkit=11.6
+```
+
+# Usage
+
+The usage is simple:
+
+```pycon
+>>> from glcm_cupy import GLCM
+>>> import numpy as np
+>>> from PIL import Image
+>>> ar = np.asarray(Image.open("image.jpg"))
+>>> ar.shape
+(1080, 1920, 3)
+>>> g = GLCM(...).run(ar)
+>>> g.shape
+(1074, 1914, 3, 8)
+```
+
+The last dimension of `g` is the GLCM Features.
+
+To retrieve a specific GLCM Feature:
+
+```pycon
+>>> from glcm_cupy import CONTRAST
+>>> g[..., CONTRAST].shape
+(1074, 1914, 3)
+```
+
+You may also consider simply `glcm` if you're not reusing `GLCM()`
+```pycon
+>>> from glcm_cupy import glcm
+>>> g = glcm(ar, ...)
+```
+
+## Features
+
+These are the features implemented.
+
+- `HOMOGENEITY = 0`
+- `CONTRAST = 1`
+- `ASM = 2`
+- `MEAN_I = 3`
+- `MEAN_J = 4`
+- `VAR_I = 5`
+- `VAR_J = 6`
+- `CORRELATION = 7`
+
+Don't see one you need? Raise an issue, I'll (hopefully) add it.
+
+## Radius & Step Size
+
+- The radius defines the window radius for each GLCM window.
+- The step size defines the distance between each window.
+  - If it's diagonal, it treats a diagonal step as 1. It's not the euclidean distance.
+
+## Binning
+
+To reduce GLCM processing time, you can specify `bin_from` & `bin_to`.
+
+This will bin the image from a range to another.
+
+I highly recommend using this to reduce time taken before raising it.
+
+E.g.
+
+> I have an RGB image with a max value of 255.
+> 
+> I limit the max value to 31. This reduces the processing time.
+> 
+> `GLCM(..., bin_from=256, bin_to=32).run(ar)`
+
+The lower the max value, the smaller the GLCM required. Thus allowing for
+more GLCMs to run concurrently.
+
+## Direction
+
+By default we have the following directions to run GLCM on.
+
+- East: `Direction.EAST`
+- South East: `Direction.SOUTH_EAST`
+- South: `Direction.SOUTH`
+- South West: `Direction.SOUTH_WEST`
+
+For each direction, the GLCM will be bi-directional.
+
+We can specify only certain directions here.
+
+```pycon
+>>> from glcm_cupy import GLCM
+>>> GLCM()
+>>> g = GLCM(directions=(Direction.SOUTH_WEST, Direction.SOUTH))
+```
+
+The result of these directions will be averaged together.
+
+# Notes
+
+> Q: Why did my image shrink?
+> 
+> The image shrunk due to `step_size` & `radius`.
+> 
+> The amount of shrink per XY Dimension is
+> `size - 2 * step_size - 2 * radius`
+
+> Q: What's the difference between this and `glcmbin5`?
+> 
+> This is the faster one, and easier to use.
+> I highly recommend avoiding `glcmbin5` as it has C++, which means you need to compile manually.
+> 
+> It's the first version of GLCM I made.
+
+## CUDA Notes
 
 ### Why is the kernel split into 4?
 
 The kernel is split into 4 sections
-1) GLCM Creation 
+
+1) GLCM Creation
 2) Features (ASM, Contrast, Homogeneity, GLCM Mean I, GLCM Mean J)
 3) Features (GLCM Variance I, GLCM Variance J)
 4) Features (GLCM Correlation)
@@ -19,23 +154,25 @@ The reason why it's split is due to (2) being reliant on (1), and (3) on (2), ..
 There are some other solutions tried
 
 1) `__syncthreads()` will not work as we require to sync all blocks.
-   1) We can't put all calculations in a block due to the thread limit of 512, 1024, 2048.
-   2) We require 256 * 256 threads minimum to support a GLCM of max value 255. 
+    1) We can't put all calculations in a block due to the thread limit of 512, 1024, 2048.
+    2) We require 256 * 256 threads minimum to support a GLCM of max value 255.
 2) **Cooperative Groups** imposes a 24 block limit.
 
 Thus, the best solution is to split the kernel.
 
 ### Atomic Add
 
-Threads cannot write to a single pointer in parallel, information will be overwritten and lost.
-This is the **Race Condition**.
+Threads cannot write to a single pointer in parallel, information will be overwritten and lost. This is the **Race
+Condition**.
 
-In order to avoid this, we use [**Atomic Functions**](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomic-functions).
+In order to avoid this, we use [**Atomic
+Functions**](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomic-functions).
 
 > ... it is guaranteed to be performed without interference from other threads
 
 ### Custom Atomic Add
 
-Currently `atomicAdd()` doesn't have the signature to support `uint8` or `unsigned char`. We get
-this implementation from this [**StackOverflow Answer**](https://stackoverflow.com/questions/5447570/cuda-atomic-operations-on-unsigned-chars)
+Currently `atomicAdd()` doesn't have the signature to support `uint8` or `unsigned char`. We get this implementation
+from this [**StackOverflow
+Answer**](https://stackoverflow.com/questions/5447570/cuda-atomic-operations-on-unsigned-chars)
 
diff --git a/definitions.py b/definitions.py
deleted file mode 100644
index 3e54954..0000000
--- a/definitions.py
+++ /dev/null
@@ -1,2 +0,0 @@
-import pathlib
-ROOT_DIR = pathlib.Path(__file__).parent.absolute()
\ No newline at end of file
diff --git a/dist/glcm_cupy-0.1.0-py3-none-any.whl b/dist/glcm_cupy-0.1.0-py3-none-any.whl
new file mode 100644
index 0000000..23f7c61
Binary files /dev/null and b/dist/glcm_cupy-0.1.0-py3-none-any.whl differ
diff --git a/dist/glcm_cupy-0.1.0.tar.gz b/dist/glcm_cupy-0.1.0.tar.gz
new file mode 100644
index 0000000..c6f85a5
Binary files /dev/null and b/dist/glcm_cupy-0.1.0.tar.gz differ
diff --git a/data/image.jpg b/examples/process_an_image/image.jpg
similarity index 100%
rename from data/image.jpg
rename to examples/process_an_image/image.jpg
diff --git a/examples/process_an_image/main.py b/examples/process_an_image/main.py
new file mode 100644
index 0000000..c635923
--- /dev/null
+++ b/examples/process_an_image/main.py
@@ -0,0 +1,33 @@
+""" Example Code """
+#%%
+import numpy as np
+from PIL import Image
+
+# Here, we load in the array
+ar = np.asarray(Image.open("image.jpg"))[::4,::4]
+
+#%%
+# We may use the class variant to run GLCM
+from glcm_cupy import GLCM
+
+g = GLCM(bin_from=256, bin_to=16).run(ar)
+
+#%%
+# Alternatively, use the function variant
+from glcm_cupy import glcm
+
+g = glcm(ar, bin_from=256, bin_to=16)
+
+#%%
+# We yield the features using constants defined in conf
+from glcm_cupy.conf import CONTRAST, CORRELATION, ASM
+
+g[..., CONTRAST]
+g[..., CORRELATION]
+g[..., ASM]
+#%%
+# Alternatively, since these constants are simply integers
+g[..., 0]
+g[..., 1]
+g[..., 2]
+
diff --git a/glcm.py b/glcm.py
deleted file mode 100644
index 2c3fbb2..0000000
--- a/glcm.py
+++ /dev/null
@@ -1,332 +0,0 @@
-from __future__ import annotations
-
-import math
-import os 
-from dataclasses import dataclass
-from typing import Tuple
-
-import cupy as cp
-import numpy as np
-from skimage.util import view_as_windows
-from tqdm import tqdm
-
-from kernel import glcm_module
-
-MAX_VALUE_SUPPORTED = 256
-NO_OF_VALUES_SUPPORTED = 256 ** 2
-MAX_RADIUS_SUPPORTED = 127
-
-MAX_THREADS = 512  # Lowest Maximum supported threads.
-
-NO_OF_FEATURES = 8
-
-# For a 10000 x 256 x 256 GLCM, you need ~ 600MB of memory.
-PARTITION_SIZE = 10000
-
-
-@dataclass
-class GLCM:
-    """
-    
-    Args:
-        max_value: Maximum value of the image, default 256
-        step_size: Step size of the window
-        radius: Radius of the windows
-        bins: Bin reduction. If None, then no reduction is done
-
-    """
-
-    step_size: int = 1
-    radius: int = 2
-    bins_from: int = 256
-    bins: int = 256
-
-    threads = MAX_VALUE_SUPPORTED + 1
-
-    HOMOGENEITY = 0
-    CONTRAST = 1
-    ASM = 2
-    MEAN_I = 3
-    MEAN_J = 4
-    VAR_I = 5
-    VAR_J = 6
-    CORRELATION = 7
-
-    @property
-    def diameter(self):
-        return self.radius * 2 + 1
-
-    @property
-    def no_of_values(self):
-        return self.diameter ** 2
-
-    def __post_init__(self):
-        self.i_flat = cp.zeros((self.diameter ** 2,), dtype=cp.uint8)
-        self.j_flat = cp.zeros((self.diameter ** 2,), dtype=cp.uint8)
-
-        if not 0 <= self.radius < MAX_RADIUS_SUPPORTED:
-            f"Radius {self.radius} should be in [0, {MAX_RADIUS_SUPPORTED})"
-        if not (2 <= self.bins <= MAX_VALUE_SUPPORTED):
-            raise ValueError(
-                f"Bins {self.bins} should be in [2, {MAX_VALUE_SUPPORTED}]. "
-            )
-        if not 1 <= self.step_size:
-            raise ValueError(
-                f"Step Size {self.step_size} should be >= 1"
-                f"If bins == 256, just use None."
-            )
-
-        self.glcm_create_kernel = \
-            glcm_module.get_function('glcmCreateKernel')
-        self.glcm_feature_kernel_0 = \
-            glcm_module.get_function('glcmFeatureKernel0')
-        self.glcm_feature_kernel_1 = \
-            glcm_module.get_function('glcmFeatureKernel1')
-        self.glcm_feature_kernel_2 = \
-            glcm_module.get_function('glcmFeatureKernel2')
-
-        os.environ['CUPY_EXPERIMENTAL_SLICE_COPY'] = '1'
-
-    @staticmethod
-    def binarize(im: np.ndarray, from_bins: int, to_bins: int) -> np.ndarray:
-        """ Binarize an image from a certain bin to another
-
-        Args:
-            im: Image as np.ndarray
-            from_bins: From the Bin of input image
-            to_bins: To the Bin of output image
-
-        Returns:
-            Binarized Image
-
-        """
-        return (im.astype(np.float32) / from_bins * to_bins).astype(np.uint8)
-
-    def from_3dimage(self,
-                     im: np.ndarray) -> np.ndarray:
-        """ Generates the GLCM from a multi band image
-
-        Args:
-            im: A 3 dim image as an ndarray
-
-        Returns:
-            The GLCM Array 4dim with shape
-                rows, cols, channel, feature
-        """
-
-        glcm_chs = []
-        for ch in range(im.shape[-1]):
-            glcm_chs.append(self.from_2dimage(im[..., ch]))
-
-        return np.stack(glcm_chs, axis=2)
-
-    def from_2dimage(self,
-                     im: np.ndarray) -> np.ndarray:
-        """ Generates the GLCM from a single band image
-
-        Notes:
-            This will actively partition the processing by blocks
-            of PARTITION_SIZE.
-            This allows for a reduction in GLCM creation.
-
-        Args:
-            im: Image in np.ndarray. Cannot be in cp.ndarray
-
-        Returns:
-            The GLCM Array 3dim with shape
-                rows, cols, feature
-        """
-
-        windows_i, windows_j = \
-            self.make_windows(im, self.diameter, self.step_size)
-
-        glcm_features = cp.zeros(
-            (windows_i.shape[0], NO_OF_FEATURES),
-            dtype=cp.float32
-        )
-
-        windows_count = windows_i.shape[0]
-        for partition in tqdm(
-            range(math.ceil(windows_count / PARTITION_SIZE))):
-            windows_part_i = windows_i[
-                             (start := partition * PARTITION_SIZE):
-                             (end := (partition + 1) * PARTITION_SIZE)
-                             ]
-            windows_part_j = windows_j[start:end]
-            glcm_features[start:start + windows_part_i.shape[0]] = \
-                self._from_windows(
-                    windows_part_i,
-                    windows_part_j
-                )
-
-        return glcm_features.reshape(
-            im.shape[0] - self.radius * 2 - self.step_size,
-            im.shape[1] - self.radius * 2 - self.step_size,
-            NO_OF_FEATURES
-        )
-
-    def _from_windows(self,
-                      i: np.ndarray,
-                      j: np.ndarray):
-        """ Generate the GLCM from the I J Window
-
-        Examples:
-
-            >>> ar_0 = np.random.randint(0, 100, 10, dtype=np.uint8)
-            >>> ar_1 = np.random.randint(0, 100, 10, dtype=np.uint8)
-            >>> g = GLCM()._from_windows(ar_0[...,np.newaxis], ar_1[...,np.newaxis])
-
-        Notes:
-            i must be the same shape as j
-
-        Args:
-            i: I Window
-            j: J Window
-
-        Returns:
-            The GLCM array, of size (8,)
-
-        """
-        assert i.ndim == 2 and j.ndim == 2, \
-            "The input dimensions must be 2. " \
-            "The 1st Dim is the partitioned windows flattened, " \
-            "The 2nd is the window cells flattened"
-
-        self.glcm = cp.zeros(
-            (i.shape[0], self.bins, self.bins),
-            dtype=cp.uint8
-        )
-        self.features = cp.zeros(
-            (i.shape[0], 8),
-            dtype=cp.float32
-        )
-        assert i.shape == j.shape, \
-            f"Shape of i {i.shape} != j {j.shape}"
-
-        i = self.binarize(i, self.bins_from, self.bins)
-        j = self.binarize(j, self.bins_from, self.bins)
-
-        no_of_windows = i.shape[0]
-
-        if i.dtype != np.uint8 or j.dtype != np.uint8:
-            raise ValueError(
-                f"Image dtype must be np.uint8,"
-                f" i: {i.dtype} j: {j.dtype}"
-            )
-        self.i_flat = cp.asarray(i)
-        self.j_flat = cp.asarray(j)
-        grid = self.calculate_grid(i.shape[0], self.bins)
-
-        self.glcm_create_kernel(
-            grid=grid,
-            block=(MAX_THREADS,),
-            args=(
-                self.i_flat,
-                self.j_flat,
-                self.bins,
-                self.no_of_values,
-                no_of_windows,
-                self.glcm,
-                self.features
-            )
-        )
-        feature_args = dict(
-            grid=grid, block=(MAX_THREADS,),
-            args=(self.glcm, self.bins, self.no_of_values, no_of_windows,
-                  self.features))
-        self.glcm_feature_kernel_0(**feature_args)
-        self.glcm_feature_kernel_1(**feature_args)
-        self.glcm_feature_kernel_2(**feature_args)
-        return self.features[:no_of_windows]
-
-    @staticmethod
-    def calculate_grid(
-        window_count: int,
-        glcm_size: int,
-        thread_per_block: int = MAX_THREADS
-    ) -> Tuple[int, int]:
-        """ Calculates the required grid size
-
-        Notes:
-            There's 2 points where the number of threads
-
-        Returns:
-            The optimal minimum grid shape for GLCM
-
-        """
-        blocks_req_glcm_calc = window_count * glcm_size * glcm_size / thread_per_block
-        # TODO: Do we need to div by thread_per_block?
-        blocks_req_glcm_populate = glcm_size * window_count
-        blocks_req = max(blocks_req_glcm_calc, blocks_req_glcm_populate)
-        return (_ := int(blocks_req ** 0.5) + 1), _
-
-    @staticmethod
-    def make_windows(im: np.ndarray,
-                     diameter: int,
-                     step_size: int) -> Tuple[np.ndarray, np.ndarray]:
-        """ From a 2D image np.ndarray, convert it into GLCM IJ windows.
-
-        Examples:
-
-            Input 4 x 4 image. Radius = 1. Step Size = 1.
-
-            +-+-+-+-+         +-----+
-            |       |         | +---+-+
-            |       |  ---->  | |   | |
-            |       |         | |   | |
-            |       |         +-+---+ |
-            +-+-+-+-+           +-----+      flat
-              4 x 4           1 x 1 x 3 x 3  ----> 1 x 9
-                              +---+   +---+
-                              flat    flat
-
-            The output will be flattened on the x,y,
-
-            Input 5 x 5 image. Radius = 1. Step Size = 1.
-
-            1-2-+-+-+-+         1-----+    2-----+
-            3 4       |         | +---+-+  | +---+-+  3-----+    4-----+
-            |         |  ---->  | |   | |  | |   | |  | +---+-+  | +---+-+
-            |         |         | |   | |  | |   | |  | |   | |  | |   | |
-            |         |         +-+---+ |  +-+---+ |  | |   | |  | |   | |
-            |         |           +-----+    +-----+  +-+---+ |  +-+---+ |
-            +-+-+-+-+-+                                 +-----+    +-----+
-              4 x 4                         2 x 2 x 3 x 3 ----> 4 x 9
-                                            +---+   +---+ flat
-                                            flat    flat
-        Args:
-            im: Input Image
-            diameter: Diameter of Window
-            step_size: Step Size between ij pairs
-
-        Returns:
-            The windows I, J suitable for GLCM.
-            The first dimension: xy flat window indexes,
-            the last dimension: xy flat indexes within each window.
-
-        """
-        # This will yield a shape (window_i, window_j, row, col)
-        # E.g. 100x100 with 5x5 window -> 96, 96, 5, 5
-
-        if im.shape[0] - step_size - diameter + 1 <= 0 or \
-            im.shape[1] - step_size - diameter + 1 <= 0:
-            raise ValueError(
-                f"Step Size & Diameter exceeds size for windowing. "
-                f"im.shape[0] {im.shape[0]} "
-                f"- step_size {step_size} "
-                f"- diameter{diameter} + 1 <= 0 or"
-                f"im.shape[1] {im.shape[1]} "
-                f"- step_size {step_size} "
-                f"- diameter{diameter} + 1 <= 0 was not satisfied."
-            )
-
-        ij = view_as_windows(im, (diameter, diameter))
-        i: np.ndarray = ij[:-step_size, :-step_size]
-        j: np.ndarray = ij[step_size:, step_size:]
-
-        i = i.reshape((-1, i.shape[-2], i.shape[-1])) \
-            .reshape((i.shape[0] * i.shape[1], -1))
-        j = j.reshape((-1, j.shape[-2], j.shape[-1])) \
-            .reshape((j.shape[0] * j.shape[1], -1))
-
-        return i, j
diff --git a/glcm_cupy/__init__.py b/glcm_cupy/__init__.py
new file mode 100644
index 0000000..5645632
--- /dev/null
+++ b/glcm_cupy/__init__.py
@@ -0,0 +1,17 @@
+from .conf import HOMOGENEITY, CONTRAST, ASM, MEAN_I, MEAN_J, VAR_I, VAR_J, CORRELATION
+from .glcm import GLCM, glcm
+from .windowing import Direction
+
+__all__ = [
+    "glcm",
+    "GLCM",
+    "HOMOGENEITY",
+    "CONTRAST",
+    "ASM",
+    "MEAN_I",
+    "MEAN_J",
+    "VAR_I",
+    "VAR_J",
+    "CORRELATION",
+    "Direction"
+]
diff --git a/glcm_cupy/conf.py b/glcm_cupy/conf.py
new file mode 100644
index 0000000..f59b92a
--- /dev/null
+++ b/glcm_cupy/conf.py
@@ -0,0 +1,22 @@
+from pathlib import Path
+
+ROOT_DIR = Path(__file__).parent.parent.absolute()
+
+MAX_VALUE_SUPPORTED = 256
+NO_OF_VALUES_SUPPORTED = 256 ** 2
+
+MAX_THREADS = 512  # Lowest Maximum supported threads.
+
+NO_OF_FEATURES = 8
+
+# For a 10000 x 256 x 256 GLCM, you need ~ 600MB of memory.
+MAX_PARTITION_SIZE = 10000
+
+HOMOGENEITY = 0
+CONTRAST = 1
+ASM = 2
+MEAN_I = 3
+MEAN_J = 4
+VAR_I = 5
+VAR_J = 6
+CORRELATION = 7
\ No newline at end of file
diff --git a/glcm_cupy/glcm.py b/glcm_cupy/glcm.py
new file mode 100644
index 0000000..dffc093
--- /dev/null
+++ b/glcm_cupy/glcm.py
@@ -0,0 +1,368 @@
+from __future__ import annotations
+
+import math
+from typing import Tuple, List, Iterable
+
+import cupy as cp
+import numpy as np
+from tqdm import tqdm
+
+from glcm_cupy.conf import *
+from glcm_cupy.kernel import glcm_module
+from glcm_cupy.windowing import make_windows, im_shape_after_glcm, Direction
+
+
+def glcm(
+    im: np.ndarray,
+    step_size: int = 1,
+    radius: int = 2,
+    bin_from: int = 256,
+    bin_to: int = 256,
+    directions: Iterable[Direction] = (Direction.EAST,
+                                       Direction.SOUTH_EAST,
+                                       Direction.SOUTH,
+                                       Direction.SOUTH_WEST)):
+    """
+    Examples:
+        To scale down the image from a 128 max value to 32, we use
+        bin_from = 128, bin_to = 32.
+
+        The range will collapse from 128 to 32.
+
+        This thus optimizes the GLCM speed.
+
+    Args:
+        im: Image to Process
+        step_size: Stride Between GLCMs
+        radius: Radius of Window
+        bin_from: Binarize from.
+        bin_to: Binarize to.
+        directions: Directions to pair the windows.
+
+    Returns:
+
+
+    """
+    return GLCM(step_size, radius, bin_from, bin_to, directions).run(im)
+
+
+class GLCM:
+    def __init__(self,
+                 step_size: int = 1,
+                 radius: int = 2,
+                 bin_from: int = 256,
+                 bin_to: int = 256,
+                 directions: Iterable[Direction] = (Direction.EAST,
+                                                    Direction.SOUTH_EAST,
+                                                    Direction.SOUTH,
+                                                    Direction.SOUTH_WEST)
+                 ):
+        """ Initialize Settings for GLCM
+
+        Examples:
+            To scale down the image from a 128 max value to 32, we use
+            bin_from = 128, bin_to = 32.
+
+            The range will collapse from 128 to 32.
+
+            This thus optimizes the GLCM speed.
+
+        Args:
+            step_size: Stride Between GLCMs
+            radius: Radius of Window
+            bin_from: Binarize from.
+            bin_to: Binarize to.
+            directions: Directions to pair the windows.
+        """
+        self.step_size = step_size
+        self.radius = radius
+        self.bin_from = bin_from
+        self.bin_to = bin_to
+        self.directions = directions
+        self.progress = None
+
+        self.i_gpu = cp.zeros((self._diameter ** 2,), dtype=cp.uint8)
+        self.j_gpu = cp.zeros((self._diameter ** 2,), dtype=cp.uint8)
+
+        if self.radius < 0:
+            raise ValueError(
+                f"Radius {radius} should be > 0)"
+            )
+        if bin_to not in range(2, MAX_VALUE_SUPPORTED + 1):
+            raise ValueError(
+                f"Target Bins {bin_to} should be "
+                f"[2, {MAX_VALUE_SUPPORTED}]. "
+            )
+        if self.step_size <= 0:
+            raise ValueError(
+                f"Step Size {step_size} should be >= 1"
+            )
+
+        self.glcm_create_kernel = \
+            glcm_module.get_function('glcmCreateKernel')
+        self.glcm_feature_kernel_0 = \
+            glcm_module.get_function('glcmFeatureKernel0')
+        self.glcm_feature_kernel_1 = \
+            glcm_module.get_function('glcmFeatureKernel1')
+        self.glcm_feature_kernel_2 = \
+            glcm_module.get_function('glcmFeatureKernel2')
+
+    @property
+    def _diameter(self):
+        return self.radius * 2 + 1
+
+    def run(self, im: np.ndarray):
+        """ Executes running GLCM. Returns the GLCM Feature array
+
+        Args:
+            im: Image to process. Can be 2D or 3D.
+
+        Returns:
+            An np.ndarray of Shape,
+             3D: (rows, cols, channels, features),
+             2D: (rows, cols, features)
+
+        """
+        shape = im_shape_after_glcm(im.shape, self.step_size, self.radius)
+        glcm_cells = np.prod(shape) * len(self.directions)
+        self.progress = tqdm(total=glcm_cells, desc="GLCM Progress",
+                             unit=" Cells",
+                             unit_scale=True)
+
+        if im.ndim == 2:
+            return self._from_2dimage(im)
+        elif im.ndim == 3:
+            return self._from_3dimage(im)
+        else:
+            raise ValueError("Only 2D or 3D images allowed.")
+
+    def _from_3dimage(self, im: np.ndarray) -> np.ndarray:
+        """ Generates the GLCM from a multi band image
+
+        Args:
+            im: A 3 dim image as an ndarray
+
+        Returns:
+            The GLCM Array 4dim with shape
+                rows, cols, channel, feature
+        """
+
+        glcm_chs = []
+        for ch in range(im.shape[-1]):
+            glcm_chs.append(self._from_2dimage(im[..., ch]))
+
+        return np.stack(glcm_chs, axis=2)
+
+    def _from_2dimage(self,
+                      im: np.ndarray) -> np.ndarray:
+        """ Generates the GLCM from a single band image
+
+        Notes:
+            Partitions the image into blocks of PARTITION_SIZE,
+            reducing allocated GLCM Memory.
+
+        Args:
+            im: Image in np.ndarray.
+
+        Returns:
+            The GLCM Array 3dim with shape
+                rows, cols, feature
+        """
+
+        if im.ndim != 2:
+            raise ValueError(f"Image must be 2 dimensional. "
+                             f"im.ndim={im.ndim}")
+
+        # Both dims are xy flattened.
+        # windows.shape == [window_ix, cell_ix]
+        dirs_ij = make_windows(im,
+                               self.radius,
+                               self.step_size,
+                               self.directions)
+        dirs_ij: List[Tuple[np.ndarray, np.ndarray]]
+        glcm_features_dirs = []
+
+        # For each direction
+        for dir, dir_ij in zip(self.directions, dirs_ij):
+
+            windows_i, windows_j = dir_ij
+
+            windows_count = windows_i.shape[0]
+            glcm_features_dir = cp.zeros(
+                (windows_count, NO_OF_FEATURES),
+                dtype=cp.float32
+            )
+
+            # If 45000 windows, 10000 partition size,
+            # We have 5 partitions (10K, 10K, 10K, 10K, 5K)
+            partition_count = math.ceil(windows_count / MAX_PARTITION_SIZE)
+
+            for partition in range(partition_count):
+                with cp.cuda.Device() as dev:
+                    # As above, we have an uneven partition
+                    # Though [1,2,3][:5] == [1,2,3]
+                    # This means windows_part_i will be correctly partitioned.
+                    part_i = windows_i[
+                             (start := partition * MAX_PARTITION_SIZE):
+                             (end := (partition + 1) * MAX_PARTITION_SIZE)
+                             ]
+                    part_j = windows_j[start:end]
+
+                    # We don't need to figure out the leftover partition size
+                    # We can simply yield the length of the leftover partition
+                    windows_part_count = part_i.shape[0]
+
+                    glcm_features_dir[start:start + windows_part_count] = \
+                        self._from_windows(
+                            part_i,
+                            part_j
+                        )
+                    dev.synchronize()
+                    self.progress.update(windows_part_count)
+
+            shape = im_shape_after_glcm(im.shape,
+                                        step_size=self.step_size,
+                                        radius=self.radius)
+            glcm_features_dirs.append(
+                glcm_features_dir
+                    .reshape(shape[0], shape[1], NO_OF_FEATURES)
+                    .get()
+            )
+
+        return np.stack(glcm_features_dirs).mean(axis=0)
+
+    def _from_windows(self,
+                      i: np.ndarray,
+                      j: np.ndarray):
+        """ Generate the GLCM from the I J Window
+
+        Examples:
+
+            >>> ar_0 = np.random.randint(0, 100, 10, dtype=np.uint8)
+            >>> ar_1 = np.random.randint(0, 100, 10, dtype=np.uint8)
+            >>> g = GLCM()._from_windows(ar_0[...,np.newaxis], ar_1[...,np.newaxis])
+
+        Notes:
+            i must be the same shape as j
+            i.shape, j.shape should be [partition_size, window_size]
+
+            For example, if you want to process 100 windows at once,
+            each window with 25 cells, you should have
+            i.shape == j.shape == (100, 25)
+
+        Args:
+            i: I Window
+            j: J Window
+
+        Returns:
+            The GLCM feature array, of size [partition_size, 8]
+
+        """
+        if i.ndim != 2 or j.ndim != 2:
+            raise ValueError(
+                f"The input dimensions must be 2. "
+                f"i.ndim=={i.ndim}, j.ndim=={j.ndim}. "
+                "The 1st Dim is the partitioned windows flattened, "
+                "The 2nd is the window cells flattened"
+            )
+
+        if i.shape != j.shape:
+            raise ValueError(f"Shape of i {i.shape} != j {j.shape}")
+
+        # partition_size != MAX_PARTITION_SIZE
+        # It may be the leftover partition.
+        partition_size = i.shape[0]
+
+        self.glcm = cp.zeros((partition_size, self.bin_to, self.bin_to),
+                             dtype=cp.uint8)
+        self.features = cp.zeros((partition_size, NO_OF_FEATURES),
+                                 dtype=cp.float32)
+
+        i = self._binner(i, self.bin_from, self.bin_to)
+        j = self._binner(j, self.bin_from, self.bin_to)
+
+        no_of_windows = i.shape[0]
+        no_of_values = self._diameter ** 2
+
+        if i.dtype != np.uint8 or j.dtype != np.uint8:
+            raise ValueError(
+                f"Image dtype must be np.uint8,"
+                f" i: {i.dtype} j: {j.dtype}"
+            )
+
+        self.i_gpu = cp.asarray(i)
+        self.j_gpu = cp.asarray(j)
+
+        self.glcm_create_kernel(
+            grid=(grid := self._calc_grid_size(no_of_windows, self.bin_to)),
+            block=(MAX_THREADS,),
+            args=(
+                self.i_gpu,
+                self.j_gpu,
+                self.bin_to,
+                no_of_values,
+                no_of_windows,
+                self.glcm,
+                self.features
+            )
+        )
+
+        feature_args = dict(
+            grid=grid, block=(MAX_THREADS,),
+            args=(self.glcm,
+                  self.bin_to,
+                  no_of_values,
+                  no_of_windows,
+                  self.features)
+        )
+
+        self.glcm_feature_kernel_0(**feature_args)
+        self.glcm_feature_kernel_1(**feature_args)
+        self.glcm_feature_kernel_2(**feature_args)
+        del self.glcm
+        return self.features[:no_of_windows]
+
+    @staticmethod
+    def _calc_grid_size(
+        window_count: int,
+        glcm_size: int,
+        thread_per_block: int = MAX_THREADS
+    ) -> Tuple[int, int]:
+        """ Calculates the required grid size
+
+        Notes:
+            There's 2 points where the number of threads
+
+        Returns:
+            The optimal minimum grid shape for GLCM
+
+        """
+
+        # TODO: Do we need to div by thread_per_block?
+        # Blocks to support features
+        blocks_req_glcm_features = \
+            window_count * glcm_size * glcm_size / thread_per_block
+
+        # Blocks to support glcm populating
+        blocks_req_glcm_populate = glcm_size * window_count
+
+        # Take the maximum
+        blocks_req = max(blocks_req_glcm_features, blocks_req_glcm_populate)
+
+        # We split it to 2 dims: A -> B x B
+        return (_ := int(blocks_req ** 0.5) + 1), _
+
+    @staticmethod
+    def _binner(im: np.ndarray, bin_from: int, bin_to: int) -> np.ndarray:
+        """ Bins an image from a certain bin to another
+
+        Args:
+            im: Image as np.ndarray
+            bin_from: From the Bin of input image
+            bin_to: To the Bin of output image
+
+        Returns:
+            Binned Image
+
+        """
+        return (im.astype(np.float32) / bin_from * bin_to).astype(np.uint8)
diff --git a/glcm_cupy/kernel.py b/glcm_cupy/kernel.py
new file mode 100644
index 0000000..881bea0
--- /dev/null
+++ b/glcm_cupy/kernel.py
@@ -0,0 +1,325 @@
+import cupy as cp
+
+glcm_module = cp.RawModule(
+    code=r"""
+   #define HOMOGENEITY 0
+#define CONTRAST 1
+#define ASM 2
+#define MEAN_I 3
+#define MEAN_J 4
+#define VAR_I 5
+#define VAR_J 6
+#define CORRELATION 7
+#define NO_OF_FEATURES 8
+
+extern "C" {
+    __device__ static inline char atomicAdd(
+        unsigned char* address,
+        unsigned char val
+        )
+    {
+        // https://stackoverflow.com/questions/5447570/cuda-atomic-operations-on-unsigned-chars
+        size_t long_address_modulo = (size_t) address & 3;
+        unsigned int* base_address = (unsigned int*) ((char*) address - long_address_modulo);
+        unsigned int long_val = (unsigned int) val << (8 * long_address_modulo);
+        unsigned int long_old = atomicAdd(base_address, long_val);
+
+        if (long_address_modulo == 3) {
+            // the first 8 bits of long_val represent the char value,
+            // hence the first 8 bits of long_old represent its previous value.
+            return (char) (long_old >> 24);
+        } else {
+            // bits that represent the char value within long_val
+            unsigned int mask = 0x000000ff << (8 * long_address_modulo);
+            unsigned int masked_old = long_old & mask;
+            // isolate the bits that represent the char value within long_old, add the long_val to that,
+            // then re-isolate by excluding bits that represent the char value
+            unsigned int overflow = (masked_old + long_val) & ~mask;
+            if (overflow) {
+                atomicSub(base_address, overflow);
+            }
+            return (char) (masked_old >> 8 * long_address_modulo);
+        }
+    }
+    __global__ void glcmCreateKernel(
+        const unsigned char* windows_i,
+        const unsigned char* windows_j,
+        const int glcmSize,
+        const int noOfValues,
+        const int noOfWindows,
+        unsigned char* g,
+        float* features)
+    {
+        /**
+        =====================================
+        Definitions
+        =====================================
+
+        windows_i, windows_j = 2D Array. Shape: (noOfValues, noOfWindows)
+
+        It should be an array of windows to be used
+
+        Take for example the following input:
+
+            +---------+  +---------+  +---------+  +---------+
+            |         |  |         |  |         |  |         |
+            |  3 x 3  |  |         |  |         |  |         |
+            |         |  |         |  |         |  |         |
+            +---------+  +---------+  +---------+  +---------+
+
+            +---------+  +---------+  +---------+  +---------+
+            |         |  |         |  |         |  |         |
+            |         |  |         |  |         |  |         |
+            |         |  |         |  |         |  |         |
+            +---------+  +---------+  +---------+  +---------+
+
+            We don't require the y dimensions. Thus, [3 x 3] x [4 x 2] = 9 x 8
+
+            windows_i will thus be 9 x 8 large. Though in CuPy, we simply flatten everything anyways.m
+            It's a simpler way to represent how it works.
+
+        glcmSize = The maximum value of the input + 1, the maximum GLCM size.
+
+        noOfValues = Number of values per window. In the above, it's 3 x 3 = 9
+
+        noOfWindows = Number of windows in total. In the above, it's 4 x 2 = 8
+
+        g = Empty initialized GLCM array. Shape of (glcmSize, glcmSize, noOfWindows)
+
+        features = Empty initialized feature array. Shape of (8, noOfWindows)
+
+        **/
+
+        /**
+        =====================================
+        Thread ID Calculation
+        =====================================
+        We are not interested in matching the block/thread dim to
+        any dims of the input. We just want to ensure that we have
+        enough threads in total.
+        **/
+
+        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
+        int tid = blockId * blockDim.x + threadIdx.x;
+
+        /**
+        =====================================
+        TID to respective windows
+        =====================================
+        TID is 1D, we need to partition them to their windows
+        Each window area is `noOfValues` == diameter ^ 2
+        We have total of `noOfWindows` number of windows
+        Thus, we simply take `wid = tid / noOfValues`
+        WID: Window ID
+
+        Note integer division
+        Example:
+            1  / 50 = 0
+            49 / 50 = 0
+            50 / 50 = 1
+            51 / 50 = 1
+
+        If we have `noOfWindows == 4`, we reject any results >= 4
+        Example:
+            1   / 50 = 0 (Accept)
+            199 / 50 = 3 (Accept)
+            200 / 50 = 4 (Reject)
+            201 / 50 = 4 (Reject)
+
+        **/
+
+        const int glcmArea = glcmSize * glcmSize;
+
+        int wid_image = tid / noOfValues;
+        const static unsigned char x = 1;
+        if (tid < noOfValues * noOfWindows)
+        {
+            unsigned char row = windows_i[tid];
+            unsigned char col = windows_j[tid];
+            // Remember that the shape of GLCM is (glcmSize, glcmSize, noOfWindows)
+            atomicAdd(&(
+                g[
+                col +
+                row * glcmSize +
+                wid_image * glcmArea
+                ]), x);
+            atomicAdd(&(
+                g[
+                row +
+                col * glcmSize +
+                wid_image * glcmArea
+                ]), x);
+        }
+    }
+
+    __global__ void glcmFeatureKernel0(
+        const unsigned char* g,
+        const int glcmSize,
+        const int noOfValues,
+        const int noOfWindows,
+        float* features)
+    {
+
+        /**
+        ===================================
+        TID to respective GLCM i, j, window
+        ===================================
+
+        To calculate the window for the tid, we simply divide by the max i & j
+
+        +---------+ +---------+ +---------+
+        |         | |         | |         |
+        |  i x j  | |         | |         |
+        |         | |         | |         |
+        +---------+ +---------+ +---------+
+        Window 0    Window 1    Window 2    ...
+
+        The area for each GLCM window is glcmSize * glcmSize
+
+        wid = tid / glcmArea
+        i = (tid % glcmArea) / glcmSize
+        j = (tid % glcmArea) % glcmSize
+
+        **/
+
+        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
+        int tid = blockId * blockDim.x + threadIdx.x;
+
+        const int glcmArea = glcmSize * glcmSize;
+        const int wid = tid / glcmArea;
+        if (wid >= noOfWindows) return;
+
+        const float i = (float)((tid % glcmArea) / glcmSize);
+        const float j = (float)((tid % glcmArea) % glcmSize);
+
+        float p = (float)(g[tid]) / noOfValues;
+
+        /**
+        =====================================
+        Feature Calculation
+        =====================================
+
+        For each feature, we require a wid * NO_OF_FEATURES offset.
+
+        8 x 1 for each GLCM
+        +----------------+ +----------------+ +----------------+
+        | HOMOGENEITY    | | HOMOGENEITY    | | HOMOGENEITY    |
+        | CONTRAST       | | CONTRAST       | | CONTRAST       |
+        | ASM            | | ASM            | | ASM            |
+        | MEAN_I         | | MEAN_I         | | MEAN_I         |
+        | MEAN_J         | | MEAN_J         | | MEAN_J         |
+        | VAR_I          | | VAR_I          | | VAR_I          |
+        | VAR_J          | | VAR_J          | | VAR_J          |
+        | CORRELATION    | | CORRELATION    | | CORRELATION    |
+        +----------------+ +----------------+ +----------------+
+        Window 0           Window 1           Window 2           ...
+        **/
+
+        __syncthreads();
+
+        atomicAdd(
+            &features[HOMOGENEITY + wid * NO_OF_FEATURES],
+            p / (1 + powf((i - j), 2.0f))
+        );
+
+        atomicAdd(
+            &features[CONTRAST + wid * NO_OF_FEATURES],
+            p * powf(i - j, 2.0f)
+        );
+
+        atomicAdd(
+            &features[ASM + wid * NO_OF_FEATURES],
+            powf(p, 2.0f)
+        );
+
+        atomicAdd(
+            &features[MEAN_I + wid * NO_OF_FEATURES],
+            p * i
+        );
+
+        atomicAdd(
+            &features[MEAN_J + wid * NO_OF_FEATURES],
+            p * j
+        );
+    }
+    __global__ void glcmFeatureKernel1(
+        const unsigned char* g,
+        const int glcmSize,
+        const int noOfValues,
+        const int noOfWindows,
+        float* features)
+    {
+        /**
+        =====================================
+        Feature Calculation
+        =====================================
+
+        See above.
+        **/
+
+        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
+        int tid = blockId * blockDim.x + threadIdx.x;
+
+        const int glcmArea = glcmSize * glcmSize;
+        const int wid = tid / glcmArea;
+        if (wid >= noOfWindows) return;
+
+        const float i = (float)((tid % glcmArea) / glcmSize);
+        const float j = (float)((tid % glcmArea) % glcmSize);
+
+        float p = (float)(g[tid]) / noOfValues;
+
+        atomicAdd(
+            &features[VAR_I + wid * NO_OF_FEATURES],
+            p * powf((i - features[MEAN_I + wid * NO_OF_FEATURES]), 2.0f)
+        );
+
+        atomicAdd(
+            &features[VAR_J + wid * NO_OF_FEATURES],
+            p * powf((j - features[MEAN_J + wid * NO_OF_FEATURES]), 2.0f)
+        );
+
+    }
+
+    __global__ void glcmFeatureKernel2(
+        const unsigned char* g,
+        const int glcmSize,
+        const int noOfValues,
+        const int noOfWindows,
+        float* features)
+    {
+        /**
+        =====================================
+        Feature Calculation
+        =====================================
+
+        See above.
+        **/
+
+        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
+        int tid = blockId * blockDim.x + threadIdx.x;
+
+        const int glcmArea = glcmSize * glcmSize;
+        const int wid = tid / glcmArea;
+        if (wid >= noOfWindows) return;
+
+        // As we invert Variance, they should never be 0.
+        if (features[VAR_I + wid * NO_OF_FEATURES] == 0 ||
+            features[VAR_J + wid * NO_OF_FEATURES] == 0) return;
+
+        const float i = (float)((tid % glcmArea) / glcmSize);
+        const float j = (float)((tid % glcmArea) % glcmSize);
+
+        float p = (float)(g[tid]) / noOfValues;
+
+        atomicAdd(
+            &features[CORRELATION + wid * NO_OF_FEATURES],
+            p * (i - features[MEAN_I + wid * NO_OF_FEATURES])
+              * (j - features[MEAN_J + wid * NO_OF_FEATURES])
+              * rsqrtf(features[VAR_I + wid * NO_OF_FEATURES]
+                     * features[VAR_J + wid * NO_OF_FEATURES])
+        );
+    }
+}
+"""
+)
diff --git a/glcm_cupy/windowing.py b/glcm_cupy/windowing.py
new file mode 100644
index 0000000..e0485e6
--- /dev/null
+++ b/glcm_cupy/windowing.py
@@ -0,0 +1,166 @@
+from __future__ import annotations
+
+from enum import Enum
+from typing import Tuple, List, Iterable
+
+import numpy as np
+from skimage.util import view_as_windows
+
+
+class Direction(Enum):
+    EAST = 0
+    SOUTH_EAST = 1
+    SOUTH = 2
+    SOUTH_WEST = 3
+
+
+def make_windows(im: np.ndarray,
+                 radius: int,
+                 step_size: int,
+                 directions: Iterable[Direction] = \
+                     (Direction.EAST,
+                      Direction.SOUTH_EAST,
+                      Direction.SOUTH,
+                      Direction.SOUTH_WEST)
+                 ) -> List[Tuple[np.ndarray, np.ndarray]]:
+    """ From a 2D image np.ndarray, convert it into GLCM IJ windows.
+
+    Examples:
+
+        Input 4 x 4 image. Radius = 1. Step Size = 1.
+
+        +-+-+-+-+         +-----+
+        |       |         | +---+-+
+        |       |  ---->  | |   | |
+        |       |         | |   | |
+        |       |         +-+---+ |
+        +-+-+-+-+           +-----+      flat
+          4 x 4           1 x 1 x 3 x 3  ----> 1 x 9
+                          +---+   +---+
+                          flat    flat
+
+        The output will be flattened on the x,y,
+
+        Input 5 x 5 image. Radius = 1. Step Size = 1.
+
+        1-2-+-+-+-+         1-----+    2-----+
+        3 4       |         | +---+-+  | +---+-+  3-----+    4-----+
+        |         |  ---->  | |   | |  | |   | |  | +---+-+  | +---+-+
+        |         |         | |   | |  | |   | |  | |   | |  | |   | |
+        |         |         +-+---+ |  +-+---+ |  | |   | |  | |   | |
+        |         |           +-----+    +-----+  +-+---+ |  +-+---+ |
+        +-+-+-+-+-+                                 +-----+    +-----+
+          4 x 4                         2 x 2 x 3 x 3 ----> 4 x 9
+                                        +---+   +---+ flat
+                                        flat    flat
+    Args:
+        im: Input Image
+        radius: Radius of Window
+        step_size: Step Size between ij pairs
+        directions: Directions to pair the windows.
+
+    Returns:
+        A List of I, J windows based on the directions.
+
+        For each window:
+            The first dimension: xy flat window indexes,
+            the last dimension: xy flat indexes within each window.
+
+    """
+
+    # This will yield a shape (window_i, window_j, row, col)
+    # E.g. 100x100 with 5x5 window -> 96, 96, 5, 5
+    shape = im_shape_after_glcm(im.shape, step_size=step_size, radius=radius)
+    if shape[0] <= 0 or shape[1] <= 0:
+        raise ValueError(
+            f"Step Size & Diameter exceeds size for windowing. "
+            f"im.shape[0] {im.shape[0]} "
+            f"- 2 * step_size {step_size} "
+            f"- 2 * radius {radius} <= 0 or"
+            f"im.shape[1] {im.shape[1]} "
+            f"- 2 * step_size {step_size} "
+            f"- 2 * radius {radius} + 1 <= 0 was not satisfied."
+        )
+
+    diameter = radius * 2 + 1
+    ij = view_as_windows(im, (diameter , diameter))
+
+    ijs: List[Tuple[np.ndarray, np.ndarray]] = []
+
+    for direction in directions:
+        i, j = pair_windows(ij, step_size=step_size, direction=direction)
+
+        i = i.reshape((-1, i.shape[-2], i.shape[-1])) \
+            .reshape((i.shape[0] * i.shape[1], -1))
+        j = j.reshape((-1, j.shape[-2], j.shape[-1])) \
+            .reshape((j.shape[0] * j.shape[1], -1))
+        ijs.append((i, j))
+
+    return ijs
+
+
+def pair_windows(ij: np.ndarray,
+                 step_size: int,
+                 direction: Direction):
+    """ Pairs the ij windows
+
+    Notes:
+
+        For an image:
+
+        +-----------+
+        |           |
+        |  +-----+  |
+        |  |     |  |
+        |  |     |  |
+        |  +-----+  |
+        |           |
+        +-----------+
+                  <-> Distance = Step Size
+
+        The i window will always be the one in the middle
+        We move j window around the outer skirts to define the pair.
+
+    Args:
+        ij: The ij output from make_windows
+        step_size: Step Size between pairs
+        direction: Direction to pair
+
+    Returns:
+
+    """
+    i = ij[step_size:-step_size, step_size:-step_size]
+    if direction == Direction.EAST:
+        j = ij[step_size:-step_size, step_size * 2:]
+    elif direction == Direction.SOUTH:
+        j = ij[step_size * 2:, step_size:-step_size]
+    elif direction == Direction.SOUTH_WEST:
+        j = ij[step_size * 2:, :-step_size * 2]
+    elif direction == Direction.SOUTH_EAST:
+        j = ij[step_size * 2:, step_size * 2:]
+    else:
+        raise ValueError("direction must be of class Direction.")
+
+    return i, j
+
+
+def im_shape_after_glcm(im_shape: Tuple,
+                        step_size: int,
+                        radius: int):
+    """ Calculate the image shape after GLCM
+
+    Args:
+        im_shape: Image Shape before GLCM
+        step_size: Step Size
+        radius: Radius of GLCM
+
+    Returns:
+        Shape of Image after GLCM
+    """
+
+    if len(im_shape) < 2:
+        raise ValueError(f"Image must be at least 2 dims.")
+
+    return (im_shape[0] - 2 * step_size - 2 * radius,
+            im_shape[1] - 2 * step_size - 2 * radius,
+            *im_shape[2:])
diff --git a/kernel.py b/kernel.py
deleted file mode 100644
index 8a3bfa6..0000000
--- a/kernel.py
+++ /dev/null
@@ -1,568 +0,0 @@
-import cupy as cp
-
-glcm_module = cp.RawModule(
-    code=r"""
-    #define HOMOGENEITY 0 
-    #define CONTRAST 1    
-    #define ASM 2         
-    #define MEAN_I 3      
-    #define MEAN_J 4      
-    #define VAR_I 5      
-    #define VAR_J 6       
-    #define CORRELATION 7 
-
-    extern "C" {
-        __device__ static inline char atomicAdd(
-            unsigned char* address,
-            unsigned char val
-            ) 
-        {
-            // https://stackoverflow.com/questions/5447570/cuda-atomic-operations-on-unsigned-chars
-            size_t long_address_modulo = (size_t) address & 3;
-            unsigned int* base_address = (unsigned int*) ((char*) address - long_address_modulo);
-            unsigned int long_val = (unsigned int) val << (8 * long_address_modulo);
-            unsigned int long_old = atomicAdd(base_address, long_val);
-
-            if (long_address_modulo == 3) {
-                // the first 8 bits of long_val represent the char value,
-                // hence the first 8 bits of long_old represent its previous value.
-                return (char) (long_old >> 24);
-            } else {
-                // bits that represent the char value within long_val
-                unsigned int mask = 0x000000ff << (8 * long_address_modulo);
-                unsigned int masked_old = long_old & mask;
-                // isolate the bits that represent the char value within long_old, add the long_val to that,
-                // then re-isolate by excluding bits that represent the char value
-                unsigned int overflow = (masked_old + long_val) & ~mask;
-                if (overflow) {
-                    atomicSub(base_address, overflow);
-                }
-                return (char) (masked_old >> 8 * long_address_modulo);
-            }
-        }
-    }
-    __global__ void glcmCreateKernel(
-        const unsigned char* windows_i,
-        const unsigned char* windows_j,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        unsigned char* g,
-        float* features)
-    {
-        /**
-        =====================================
-        Definitions
-        =====================================
-
-        windows_i, windows_j = 2D Array. Shape: (noOfValues, noOfWindows)
-
-        It should be an array of windows to be used
-
-        Take for example the following input:
-
-            +---------+  +---------+  +---------+  +---------+
-            |         |  |         |  |         |  |         |
-            |  3 x 3  |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            +---------+  +---------+  +---------+  +---------+
-
-            +---------+  +---------+  +---------+  +---------+
-            |         |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            +---------+  +---------+  +---------+  +---------+
-
-            We don't require the y dimensions. Thus, [3 x 3] x [4 x 2] = 9 x 8
-
-            windows_i will thus be 9 x 8 large. Though in CuPy, we simply flatten everything anyways.m
-            It's a simpler way to represent how it works.
-
-        glcmSize = The maximum value of the input + 1, the maximum GLCM size.
-
-        noOfValues = Number of values per window. In the above, it's 3 x 3 = 9
-
-        noOfWindows = Number of windows in total. In the above, it's 4 x 2 = 8
-
-        g = Empty initialized GLCM array. Shape of (glcmSize, glcmSize, noOfWindows)
-
-        features = Empty initialized feature array. Shape of (8, noOfWindows)
-
-        **/
-
-        /**
-        =====================================
-        Thread ID Calculation
-        =====================================
-        We are not interested in matching the block/thread dim to
-        any dims of the input. We just want to ensure that we have
-        enough threads in total.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        /**
-        =====================================
-        TID to respective windows
-        =====================================
-        TID is 1D, we need to partition them to their windows
-        Each window area is `noOfValues` == diameter ^ 2
-        We have total of `noOfWindows` number of windows
-        Thus, we simply take `wid = tid / noOfValues`
-        WID: Window ID
-
-        Note integer division
-        Example:
-            1  / 50 = 0
-            49 / 50 = 0
-            50 / 50 = 1
-            51 / 50 = 1
-
-        If we have `noOfWindows == 4`, we reject any results >= 4
-        Example:
-            1   / 50 = 0 (Accept)
-            199 / 50 = 3 (Accept)
-            200 / 50 = 4 (Reject)
-            201 / 50 = 4 (Reject)
-
-        **/
-
-        const int glcmArea = glcmSize * glcmSize;
-
-        int wid_image = tid / noOfValues;
-        const static unsigned char x = 1;
-        if (tid < noOfValues * noOfWindows)
-        {
-            int tid = threadIdx.x + blockIdx.x * blockDim.x;
-            unsigned char x = 1;
-            // Prevent OOB
-            if (tid < noOfValues){
-                unsigned char row = window_i[tid];
-                unsigned char col = window_j[tid];
-                atomicAdd(&(g[col + row * (maxValue + 1)]), x);
-            }
-            return (char) (masked_old >> 8 * long_address_modulo);
-        }
-    }
-
-    __global__ void glcmFeatureKernel0(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-
-        /**
-        ===================================
-        TID to respective GLCM i, j, window
-        ===================================
-
-        To calculate the window for the tid, we simply divide by the max i & j
-
-        +---------+ +---------+ +---------+
-        |         | |         | |         |
-        |  i x j  | |         | |         |
-        |         | |         | |         |
-        +---------+ +---------+ +---------+
-        Window 0    Window 1    Window 2    ...
-
-        The area for each GLCM window is glcmSize * glcmSize
-
-        wid = tid / glcmArea
-        i = (tid % glcmArea) / glcmSize
-        j = (tid % glcmArea) % glcmSize
-
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        float p = (float)(g[tid]) / noOfValues;
-
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        For each feature, we require a wid * NO_OF_FEATURES offset.
-
-        8 x 1 for each GLCM
-        +----------------+ +----------------+ +----------------+
-        | HOMOGENEITY    | | HOMOGENEITY    | | HOMOGENEITY    |
-        | CONTRAST       | | CONTRAST       | | CONTRAST       |
-        | ASM            | | ASM            | | ASM            |
-        | MEAN_I         | | MEAN_I         | | MEAN_I         |
-        | MEAN_J         | | MEAN_J         | | MEAN_J         |
-        | VAR_I          | | VAR_I          | | VAR_I          |
-        | VAR_J          | | VAR_J          | | VAR_J          |
-        | CORRELATION    | | CORRELATION    | | CORRELATION    |
-        +----------------+ +----------------+ +----------------+
-        Window 0           Window 1           Window 2           ...
-        **/
-
-        __syncthreads();
-
-        atomicAdd(
-            &features[HOMOGENEITY + wid * NO_OF_FEATURES],
-            p / (1 + powf((i - j), 2.0f))
-        );
-
-        atomicAdd(
-            &features[CONTRAST + wid * NO_OF_FEATURES],
-            p * powf(i - j, 2.0f)
-        );
-
-        atomicAdd(
-            &features[ASM + wid * NO_OF_FEATURES],
-            powf(p, 2.0f)
-        );
-
-        atomicAdd(
-            &features[MEAN_I + wid * NO_OF_FEATURES],
-            p * i
-        );
-
-        atomicAdd(
-            &features[MEAN_J + wid * NO_OF_FEATURES],
-            p * j
-        );
-    }
-    __global__ void glcmFeatureKernel1(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        See above.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        float p = (float)(g[tid]) / noOfValues;
-
-        atomicAdd(
-            &features[VAR_I + wid * NO_OF_FEATURES],
-            p * powf((i - features[MEAN_I + wid * NO_OF_FEATURES]), 2.0f)
-        );
-
-        atomicAdd(
-            &features[VAR_J + wid * NO_OF_FEATURES],
-            p * powf((j - features[MEAN_J + wid * NO_OF_FEATURES]), 2.0f)
-        );
-
-    }
-
-    __global__ void glcmFeatureKernel2(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        See above.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        // As we invert Variance, they should never be 0.
-        if (features[VAR_I + wid * NO_OF_FEATURES] == 0 ||
-            features[VAR_J + wid * NO_OF_FEATURES] == 0) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        windows_i, windows_j = 2D Array. Shape: (noOfValues, noOfWindows)
-
-        It should be an array of windows to be used
-
-        Take for example the following input:
-
-            +---------+  +---------+  +---------+  +---------+
-            |         |  |         |  |         |  |         |
-            |  3 x 3  |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            +---------+  +---------+  +---------+  +---------+
-
-            +---------+  +---------+  +---------+  +---------+
-            |         |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            |         |  |         |  |         |  |         |
-            +---------+  +---------+  +---------+  +---------+
-
-            We don't require the y dimensions. Thus, [3 x 3] x [4 x 2] = 9 x 8
-
-            windows_i will thus be 9 x 8 large. Though in CuPy, we simply flatten everything anyways.m
-            It's a simpler way to represent how it works.
-
-        glcmSize = The maximum value of the input + 1, the maximum GLCM size.
-
-        noOfValues = Number of values per window. In the above, it's 3 x 3 = 9
-
-        noOfWindows = Number of windows in total. In the above, it's 4 x 2 = 8
-
-        g = Empty initialized GLCM array. Shape of (glcmSize, glcmSize, noOfWindows)
-
-        features = Empty initialized feature array. Shape of (8, noOfWindows)
-
-        **/
-
-        /**
-        =====================================
-        Thread ID Calculation
-        =====================================
-        We are not interested in matching the block/thread dim to
-        any dims of the input. We just want to ensure that we have
-        enough threads in total.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        /**
-        =====================================
-        TID to respective windows
-        =====================================
-        TID is 1D, we need to partition them to their windows
-        Each window area is `noOfValues` == diameter ^ 2
-        We have total of `noOfWindows` number of windows
-        Thus, we simply take `wid = tid / noOfValues`
-        WID: Window ID
-
-        Note integer division
-        Example:
-            1  / 50 = 0
-            49 / 50 = 0
-            50 / 50 = 1
-            51 / 50 = 1
-
-        If we have `noOfWindows == 4`, we reject any results >= 4
-        Example:
-            1   / 50 = 0 (Accept)
-            199 / 50 = 3 (Accept)
-            200 / 50 = 4 (Reject)
-            201 / 50 = 4 (Reject)
-
-        **/
-
-        const int glcmArea = glcmSize * glcmSize;
-
-        int wid_image = tid / noOfValues;
-        const static unsigned char x = 1;
-        if (tid < noOfValues * noOfWindows)
-        {
-            unsigned char row = windows_i[tid];
-            unsigned char col = windows_j[tid];
-            // Remember that the shape of GLCM is (glcmSize, glcmSize, noOfWindows)
-            atomicAdd(&(
-                g[
-                col +
-                row * glcmSize +
-                wid_image * glcmArea
-                ]), x);
-        }
-    }
-
-    __global__ void glcm_1(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-
-        /**
-        ===================================
-        TID to respective GLCM i, j, window
-        ===================================
-
-        To calculate the window for the tid, we simply divide by the max i & j
-
-        +---------+ +---------+ +---------+
-        |         | |         | |         |
-        |  i x j  | |         | |         |
-        |         | |         | |         |
-        +---------+ +---------+ +---------+
-        Window 0    Window 1    Window 2    ...
-
-        The area for each GLCM window is glcmSize * glcmSize
-
-        wid = tid / glcmArea
-        i = (tid % glcmArea) / glcmSize
-        j = (tid % glcmArea) % glcmSize
-
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        float p = (float)(g[tid]) / noOfValues;
-
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        For each feature, we require a wid * NO_OF_FEATURES offset.
-
-        8 x 1 for each GLCM
-        +----------------+ +----------------+ +----------------+
-        | HOMOGENEITY    | | HOMOGENEITY    | | HOMOGENEITY    |
-        | CONTRAST       | | CONTRAST       | | CONTRAST       |
-        | ASM            | | ASM            | | ASM            |
-        | MEAN_I         | | MEAN_I         | | MEAN_I         |
-        | MEAN_J         | | MEAN_J         | | MEAN_J         |
-        | VAR_I          | | VAR_I          | | VAR_I          |
-        | VAR_J          | | VAR_J          | | VAR_J          |
-        | CORRELATION    | | CORRELATION    | | CORRELATION    |
-        +----------------+ +----------------+ +----------------+
-        Window 0           Window 1           Window 2           ...
-        **/
-
-        __syncthreads();
-
-        atomicAdd(
-            &features[HOMOGENEITY + wid * NO_OF_FEATURES],
-            p / (1 + powf((i - j), 2.0f))
-        );
-
-        atomicAdd(
-            &features[CONTRAST + wid * NO_OF_FEATURES],
-            p * powf(i - j, 2.0f)
-        );
-
-        atomicAdd(
-            &features[ASM + wid * NO_OF_FEATURES],
-            powf(p, 2.0f)
-        );
-
-        atomicAdd(
-            &features[MEAN_I + wid * NO_OF_FEATURES],
-            p * i
-        );
-
-        atomicAdd(
-            &features[MEAN_J + wid * NO_OF_FEATURES],
-            p * j
-        );
-    }
-    __global__ void glcm_2(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        See above.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        float p = (float)(g[tid]) / noOfValues;
-
-        atomicAdd(
-            &features[VAR_I + wid * NO_OF_FEATURES],
-            p * powf((i - features[MEAN_I + wid * NO_OF_FEATURES]), 2.0f)
-        );
-
-        atomicAdd(
-            &features[VAR_J + wid * NO_OF_FEATURES],
-            p * powf((j - features[MEAN_J + wid * NO_OF_FEATURES]), 2.0f)
-        );
-
-    }
-
-    __global__ void glcm_3(
-        const unsigned char* g,
-        const int glcmSize,
-        const int noOfValues,
-        const int noOfWindows,
-        float* features)
-    {
-        /**
-        =====================================
-        Feature Calculation
-        =====================================
-
-        See above.
-        **/
-
-        int blockId = blockIdx.y * gridDim.x + blockIdx.x;
-        int tid = blockId * blockDim.x + threadIdx.x;
-
-        const int glcmArea = glcmSize * glcmSize;
-        const int wid = tid / glcmArea;
-        if (wid >= noOfWindows) return;
-
-        // As we invert Variance, they should never be 0.
-        if (features[VAR_I + wid * NO_OF_FEATURES] == 0 ||
-            features[VAR_J + wid * NO_OF_FEATURES] == 0) return;
-
-        const float i = (float)((tid % glcmArea) / glcmSize);
-        const float j = (float)((tid % glcmArea) % glcmSize);
-
-        float p = (float)(g[tid]) / noOfValues;
-
-        atomicAdd(
-            &features[CORRELATION + wid * NO_OF_FEATURES],
-            p * (i - features[MEAN_I + wid * NO_OF_FEATURES])
-              * (j - features[MEAN_J + wid * NO_OF_FEATURES])
-              * rsqrtf(features[VAR_I + wid * NO_OF_FEATURES]
-                     * features[VAR_J + wid * NO_OF_FEATURES])
-        );
-    }
-}
-"""
-)
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..ed06a4f
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,22 @@
+[tool.poetry]
+name = "glcm_cupy"
+version = "0.1.0"
+description = "Binned GLCM 5 Features implemented in CuPy"
+authors = ["Evening"]
+license = "MIT"
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Environment :: GPU :: NVIDIA CUDA",
+    "Programming Language :: Python :: 3.8"
+]
+
+[tool.poetry.dependencies]
+python = ">=3.8"
+scikit-image = "~=0.18.3"
+
+[tool.poetry.dev-dependencies]
+scikit-image = "~=0.18.3"
+
+[build-system]
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
diff --git a/requirements.txt b/requirements.txt
index 42304d1..0684498 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,3 @@
 cupy~=10.2.0
 scikit-image~=0.18.3
+tqdm
\ No newline at end of file
diff --git a/kernel.cu b/rsc/kernel.cu
similarity index 97%
rename from kernel.cu
rename to rsc/kernel.cu
index cbde429..4b8472e 100644
--- a/kernel.cu
+++ b/rsc/kernel.cu
@@ -37,7 +37,7 @@ extern "C" {
             return (char) (masked_old >> 8 * long_address_modulo);
         }
     }
-    __global__ void glcm_0(
+    __global__ void glcmCreateKernel(
         const unsigned char* windows_i,
         const unsigned char* windows_j,
         const int glcmSize,
@@ -103,7 +103,7 @@ extern "C" {
         TID to respective windows
         =====================================
         TID is 1D, we need to partition them to their windows
-        Each window area is `noOfValues` == diameter ^ 2
+        Each window area is `noOfValues` == _diameter ^ 2
         We have total of `noOfWindows` number of windows
         Thus, we simply take `wid = tid / noOfValues`
         WID: Window ID
@@ -142,7 +142,7 @@ extern "C" {
         }
     }
 
-    __global__ void glcm_1(
+    __global__ void glcmFeatureKernel0(
         const unsigned char* g,
         const int glcmSize,
         const int noOfValues,
@@ -232,7 +232,7 @@ extern "C" {
             p * j
         );
     }
-    __global__ void glcm_2(
+    __global__ void glcmFeatureKernel1(
         const unsigned char* g,
         const int glcmSize,
         const int noOfValues,
@@ -271,7 +271,7 @@ extern "C" {
 
     }
 
-    __global__ void glcm_3(
+    __global__ void glcmFeatureKernel2(
         const unsigned char* g,
         const int glcmSize,
         const int noOfValues,
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..651be72
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,10 @@
+[metadata]
+long_description = file: README.md
+long_description_content_type = text/markdown
+url = https://github.com/eve-ning
+;project_urls =
+;    Bug Tracker = https://github.com/pypa/sampleproject/issues
+classifiers =
+    Programming Language :: Python :: 3
+    License :: OSI Approved :: MIT License
+    Operating System :: OS Independent


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
deleted file mode 100644
index 18a2a0c..0000000
--- a/tests/conftest.py
+++ /dev/null
@@ -1,8 +0,0 @@
-import os
-
-import pytest
-
-
-@pytest.fixture(autouse=True)
-def set_env():
-    os.environ['CUPY_EXPERIMENTAL_SLICE_COPY'] = '1'
diff --git a/tests/integration_tests/conftest.py b/tests/integration_tests/conftest.py
index 64bb3ce..2276939 100644
--- a/tests/integration_tests/conftest.py
+++ b/tests/integration_tests/conftest.py
@@ -2,7 +2,7 @@
 import pytest
 from PIL import Image
 
-from definitions import ROOT_DIR
+from glcm_cupy.conf import ROOT_DIR
 
 TEST_SIZE = 25
 
diff --git a/tests/integration_tests/test_glcm.py b/tests/integration_tests/test_glcm.py
index e0e7f4a..6f02165 100644
--- a/tests/integration_tests/test_glcm.py
+++ b/tests/integration_tests/test_glcm.py
@@ -1,16 +1,16 @@
 import numpy as np
 
-from glcm import GLCM
+from glcm_cupy import GLCM
 
 
 def test_from_3dimage(ar_3d):
     """ Tests using a 3D Image """
-    GLCM().from_3dimage(ar_3d)
+    GLCM()._from_3dimage(ar_3d)
 
 
 def test_from_2dimage(ar_2d):
     """ Tests with a 2D Image (1 Channel) """
-    GLCM().from_2dimage(ar_2d)
+    GLCM()._from_2dimage(ar_2d)
 
 
 def test__from_windows():
diff --git a/tests/integration_tests/test_image.py b/tests/integration_tests/test_image.py
index a1c9657..07d16c3 100644
--- a/tests/integration_tests/test_image.py
+++ b/tests/integration_tests/test_image.py
@@ -1,11 +1,11 @@
 import numpy as np
 from PIL import Image
 
-from definitions import ROOT_DIR
-from glcm import GLCM
+from glcm_cupy.conf import ROOT_DIR
+from glcm_cupy import GLCM
 
 
 def test_image():
     img = Image.open(f"{ROOT_DIR}/data/image.jpg")
-    ar = np.asarray(img)[::5]
-    _ = GLCM().from_3dimage(ar)
+    ar = np.asarray(img)[::5, ::5]
+    _ = GLCM(bin_to=16)._from_3dimage(ar)
diff --git a/tests/unit_tests/test_glcm.py b/tests/unit_tests/test_glcm.py
index 6c5a5df..c19953d 100644
--- a/tests/unit_tests/test_glcm.py
+++ b/tests/unit_tests/test_glcm.py
@@ -5,7 +5,8 @@
 import pytest
 from pytest_mock import MockerFixture
 
-from glcm import GLCM, PARTITION_SIZE
+from glcm_cupy import GLCM
+from glcm_cupy.conf import *
 from tests.unit_tests import glcm_py
 
 
@@ -37,14 +38,14 @@ def test_glcm_from_windows(i, j):
 
     # The sum of the values, since tiled, will be scaled by no of windows.
     actual = dict(
-        homogeneity=float(g[..., GLCM.HOMOGENEITY].sum() / windows),
-        contrast=float(g[..., GLCM.CONTRAST].sum() / windows),
-        asm=float(g[..., GLCM.ASM].sum() / windows),
-        mean_i=float(g[..., GLCM.MEAN_I].sum() / windows),
-        mean_j=float(g[..., GLCM.MEAN_J].sum() / windows),
-        var_i=float(g[..., GLCM.VAR_I].sum() / windows),
-        var_j=float(g[..., GLCM.VAR_J].sum() / windows),
-        correlation=float(g[..., GLCM.CORRELATION].sum() / windows)
+        homogeneity=float(g[..., HOMOGENEITY].sum() / windows),
+        contrast=float(g[..., CONTRAST].sum() / windows),
+        asm=float(g[..., ASM].sum() / windows),
+        mean_i=float(g[..., MEAN_I].sum() / windows),
+        mean_j=float(g[..., MEAN_J].sum() / windows),
+        var_i=float(g[..., VAR_I].sum() / windows),
+        var_j=float(g[..., VAR_J].sum() / windows),
+        correlation=float(g[..., CORRELATION].sum() / windows)
     )
 
     expected = glcm_py(i, j)
@@ -91,20 +92,20 @@ def test_glcm_make_windows(
     if y_windows <= 0 or x_windows <= 0:
         # If the make windows is invalid, we assert that it throws an error
         with pytest.raises(ValueError):
-            GLCM.make_windows(im, radius * 2 + 1, step_size)
+            GLCM._make_windows(im, radius * 2 + 1, step_size)
     else:
         # Else, we assert the correct shape returns
-        windows = GLCM.make_windows(im, radius * 2 + 1, step_size)
+        windows = GLCM._make_windows(im, radius * 2 + 1, step_size)
         assert (x_windows * y_windows, x_cells * y_cells) == windows[0].shape
 
 
 @pytest.mark.parametrize(
     "windows",
     [1,
-     PARTITION_SIZE,
-     PARTITION_SIZE + 1,
-     PARTITION_SIZE * 2,
-     PARTITION_SIZE * 2 + 1]
+     MAX_PARTITION_SIZE,
+     MAX_PARTITION_SIZE + 1,
+     MAX_PARTITION_SIZE * 2,
+     MAX_PARTITION_SIZE * 2 + 1]
 )
 def test_glcm_partition(
     windows: int,
@@ -133,7 +134,7 @@ def __setitem__(self, key, value): ...
 
     mocker.patch('cupy.zeros', return_value=MockFeatures())
     mocker.patch.object(
-        glcm_instance, 'make_windows',
+        glcm_instance, '_make_windows',
         return_value=(
             np.zeros((windows, WINDOW_SIZE), dtype=np.uint8),
             np.zeros((windows, WINDOW_SIZE), dtype=np.uint8),
@@ -141,7 +142,7 @@ def __setitem__(self, key, value): ...
     )
     mock_from_windows = \
         mocker.patch.object(glcm_instance, '_from_windows', return_value=[])
-    glcm_instance.from_2dimage(np.asarray([[]]))
+    glcm_instance._from_2dimage(np.asarray([[]]))
 
     # The LHS extracts the array size calls to _from_windows
     #  We expect the windows to be split to chunks of PARTITION_SIZE
@@ -151,8 +152,8 @@ def __setitem__(self, key, value): ...
     #  The second is the remainder, however 0 will produce a [] as expected.
     assert [call[0][0].shape[0] for call in
             mock_from_windows.call_args_list] == \
-           [PARTITION_SIZE for _ in range(windows // PARTITION_SIZE)] + \
-           ([windows % PARTITION_SIZE] if windows % PARTITION_SIZE else [])
+           [MAX_PARTITION_SIZE for _ in range(windows // MAX_PARTITION_SIZE)] + \
+           ([windows % MAX_PARTITION_SIZE] if windows % MAX_PARTITION_SIZE else [])
 
 
 @pytest.mark.parametrize(
@@ -173,7 +174,7 @@ def test_glcm_binarize(bins):
     Args:
         bins: The result bins
     """
-    g = GLCM.binarize(np.asarray([0, 1, 2], dtype=np.uint8), 3, bins)
+    g = GLCM._binner(np.asarray([0, 1, 2], dtype=np.uint8), 3, bins)
 
     assert g[0] == 0 // 3
     assert g[1] == bins // 3