diff --git a/python/cubes.py b/python/cubes.py
index 1b5647a..df62b38 100644
--- a/python/cubes.py
+++ b/python/cubes.py
@@ -5,7 +5,7 @@
 from libraries.resizing import expand_cube
 from libraries.packing import pack, unpack
 from libraries.renderer import render_shapes
-from libraries.rotation import all_rotations
+from libraries.rotation import all_rotations, one_diff_rot, all_diff_rot
 
 
 def log_if_needed(n, total_n):
@@ -87,11 +87,32 @@ def get_canonical_packing(polycube: np.ndarray,
 
     """
     max_id = b'\x00'
-    for cube_rotation in all_rotations(polycube):
+    shape = polycube.shape
+
+    if shape[0] == shape[1] == shape[2]:
+        for cube_rotation in all_rotations(polycube):
+            this_id = pack(cube_rotation)
+            if this_id in known_ids:
+                return this_id
+            if this_id > max_id:
+                max_id = this_id
+        return max_id
+
+    for idx in range(0, 3):
+        if shape[idx] == shape[(idx + 1) % 3]:
+            for cube_rotation in one_diff_rot(polycube, (idx + 2) % 3, [idx, (idx + 1) % 3]):
+                this_id = pack(cube_rotation)
+                if this_id in known_ids:
+                    return this_id
+                if this_id > max_id:
+                    max_id = this_id
+            return max_id
+
+    for cube_rotation in all_diff_rot(polycube):
         this_id = pack(cube_rotation)
-        if (this_id in known_ids):
+        if this_id in known_ids:
             return this_id
-        if (this_id > max_id):
+        if this_id > max_id:
             max_id = this_id
     return max_id
 
diff --git a/python/libraries/resizing.py b/python/libraries/resizing.py
index 517f5f4..d152b0b 100644
--- a/python/libraries/resizing.py
+++ b/python/libraries/resizing.py
@@ -23,6 +23,52 @@ def crop_cube(cube: np.ndarray) -> np.ndarray:
     return cube
 
 
+def fix_axis(cube):
+    """
+    Rotate cube to make boundary sizes in descending order.
+
+    Ex : if input cube.shape = (4, 1, 2) this method rotate it to be (4, 2, 1)
+
+    Parameters:
+    cube (np.array): 3D Numpy byte array where 1 values indicate polycube positions
+
+    Returns:
+    np.array: Rotated 3D Numpy byte array equivalent to cube, but with descending size orientation
+
+    """
+    if cube.shape == tuple(sorted(cube.shape, reverse=True)):
+        return cube
+
+    if cube.shape == tuple(sorted(cube.shape)):
+        return np.rot90(cube, 1, (0, 2))
+
+    if cube.shape[0] < cube.shape[2] < cube.shape[1]:
+        return np.rot90(np.rot90(cube, 1, (0, 2)), 1, (1, 2))
+
+    if cube.shape[0] == cube.shape[2]:
+        if cube.shape[0] < cube.shape[1]:
+            return np.rot90(cube, 1, (0, 1))
+        else:
+            return np.rot90(cube, 1, (1, 2))
+
+    if cube.shape[0] == cube.shape[1] or cube.shape[1] == cube.shape[2]:
+        if cube.shape[0] < cube.shape[2]:
+            return cube
+        else:
+            return np.rot90(cube, 1, (0, 2))
+
+    if cube.shape[2] < cube.shape[0] < cube.shape[1]:
+        return np.rot90(cube, 1, (0, 1))
+
+    if cube.shape[1] < cube.shape[0] < cube.shape[2]:
+        return np.rot90(np.rot90(cube, 1, (0, 1)), 1, (1, 2))
+
+    if cube.shape[1] < cube.shape[2] < cube.shape[0]:
+        return np.rot90(cube, 1, (2, 1))
+
+    print("error", cube.shape)
+    exit(2)
+
 def expand_cube(cube: np.ndarray) -> Generator[np.ndarray, None, None]:
     """
     Expands a polycube by adding single blocks at all valid locations.
@@ -53,4 +99,4 @@ def expand_cube(cube: np.ndarray) -> Generator[np.ndarray, None, None]:
     for (x, y, z) in zip(exp[0], exp[1], exp[2]):
         new_cube = np.array(cube)
         new_cube[x, y, z] = 1
-        yield crop_cube(new_cube)
+        yield fix_axis(crop_cube(new_cube))
diff --git a/python/libraries/rotation.py b/python/libraries/rotation.py
index e20edb6..51aa5c2 100644
--- a/python/libraries/rotation.py
+++ b/python/libraries/rotation.py
@@ -4,7 +4,7 @@
 
 def all_rotations(polycube: np.ndarray) -> Generator[np.ndarray, None, None]:
     """
-    Calculates all rotations of a polycube.
+    Calculates rotations of a polycube when bounding box size in all axes are equal(Ex (3, 3, 3)).
 
     Adapted from https://stackoverflow.com/questions/33190042/how-to-calculate-all-24-rotations-of-3d-array.
     This function computes all 24 rotations around each of the axis x,y,z. It uses numpy operations to do this, to avoid unecessary copies.
@@ -36,3 +36,47 @@ def single_axis_rotation(polycube, axes):
     # rotate about axis 2, 8 rotations about axis 1
     yield from single_axis_rotation(np.rot90(polycube, axes=(0, 1)), (0, 2))
     yield from single_axis_rotation(np.rot90(polycube, -1, axes=(0, 1)), (0, 2))
+
+
+def one_diff_rot(cube, diff_axis, equal_axes):
+    """
+    Calculates rotations of a polycube when bounding box size in two axes are equal, but one is different. (Ex (2, 2, 3)).
+    Only 8 rotations can be done without breaking the bounding box.
+
+    Parameters:
+    polycube (np.array): 3D Numpy byte array where 1 values indicate polycube positions
+    diff_axis (integer): Axis which has different size
+    equal_axes (integer array): Axes which have same size
+
+    Returns:
+    generator(np.array): Yields new rotations of this cube about all axes
+
+    """
+    for _ in range(0, 4):
+        yield cube
+        cube = np.rot90(cube, 1, equal_axes)
+
+    cube = np.rot90(cube, 2, (diff_axis, equal_axes[0]))
+
+    for _ in range(0, 4):
+        yield cube
+        cube = np.rot90(cube, 1, equal_axes)
+
+
+def all_diff_rot(cube):
+    """
+    Calculates rotations of a polycube when bounding box size in all axes are different. (Ex (1, 4, 6)).
+    Only 4 rotations can be done without breaking the bounding box.
+
+    Parameters:
+    polycube (np.array): 3D Numpy byte array where 1 values indicate polycube positions
+
+    Returns:
+    generator(np.array): Yields new rotations of this cube about all axes
+
+    """
+    yield cube
+    yield np.rot90(cube, 2, (0, 1))
+    yield np.rot90(cube, 2, (0, 2))
+    yield np.rot90(cube, 2, (1, 2))
+