diff --git a/foa3d/input.py b/foa3d/input.py
index 8af4622..78c1e8e 100644
--- a/foa3d/input.py
+++ b/foa3d/input.py
@@ -84,11 +84,13 @@ def get_cli_parser():
     cli_parser.add_argument('--z-max', type=float, default=None, help='forced maximum output z-depth [μm]')
     cli_parser.add_argument('--hsv', action='store_true', default=False,
                             help='toggle HSV colormap for 3D fiber orientations')
-    cli_parser.add_argument('-o', '--odf-res', nargs='+', type=float, help='side of the fiber ODF super-voxels: '
-                                                                           'do not generate ODFs if None [μm]')
+    cli_parser.add_argument('--odf-res', nargs='+', type=float, help='side of the fiber ODF super-voxels: '
+                                                                     'do not generate ODFs if None [μm]')
     cli_parser.add_argument('--odf-deg', type=int, default=6,
                             help='degrees of the spherical harmonics series expansion (even number between 2 and 10)')
-
+    cli_parser.add_argument('-o', '--out', type=str, default=None,
+                            help='output directory')
+    
     # parse arguments
     cli_args = cli_parser.parse_args()
 
diff --git a/foa3d/output.py b/foa3d/output.py
index 345b470..c60389b 100644
--- a/foa3d/output.py
+++ b/foa3d/output.py
@@ -1,5 +1,5 @@
 from datetime import datetime
-from os import mkdir, path
+from os import makedirs, path
 
 import nibabel as nib
 import numpy as np
@@ -34,32 +34,35 @@ def create_save_dirs(img_path, img_name, cli_args, is_fiber=False):
     # get current time
     time_stamp = datetime.now().strftime("%Y%m%d-%H%M%S")
 
-    # get base path
-    base_path = path.dirname(img_path)
+    # get output path
+    out_path = cli_args.out
+    if out_path is None:
+        out_path = path.dirname(img_path)
 
     # create saving directory
-    base_dir = path.join(base_path, time_stamp + '_' + img_name)
-    save_dir = list()
-    if not path.isdir(base_dir):
-        mkdir(base_dir)
+    base_out_dir = path.join(out_path, time_stamp + '_' + img_name)
+    save_dir_lst = list()
+    if not path.isdir(base_out_dir):
+        print("Creating directory")
+        makedirs(base_out_dir)
 
     # create Frangi filter output subdirectory
     if not is_fiber:
-        frangi_dir = path.join(base_dir, 'frangi')
-        mkdir(frangi_dir)
-        save_dir.append(frangi_dir)
+        frangi_dir = path.join(base_out_dir, 'frangi')
+        makedirs(frangi_dir)
+        save_dir_lst.append(frangi_dir)
     else:
-        save_dir.append(None)
+        save_dir_lst.append(None)
 
     # create ODF analysis output subdirectory
     if cli_args.odf_res is not None:
-        odf_dir = path.join(base_dir, 'odf')
-        mkdir(odf_dir)
-        save_dir.append(odf_dir)
+        odf_dir = path.join(base_out_dir, 'odf')
+        makedirs(odf_dir)
+        save_dir_lst.append(odf_dir)
     else:
-        save_dir.append(None)
+        save_dir_lst.append(None)
 
-    return save_dir
+    return save_dir_lst
 
 
 def save_array(fname, save_dir, nd_array, px_sz=None, fmt='tiff', odi=False):
diff --git a/foa3d/utils.py b/foa3d/utils.py
index 7182e7e..c462035 100644
--- a/foa3d/utils.py
+++ b/foa3d/utils.py
@@ -468,7 +468,7 @@ def hsv_orient_cmap(vec_img):
 
     # compute the in-plane versor length
     vxy_abs = np.sqrt(np.square(vx) + np.square(vy))
-    vxy_abs = np.divide(vxy_abs, np.max(vxy_abs))
+    vxy_abs = divide_nonzero(vxy_abs, np.max(vxy_abs))
 
     # compute the in-plane angular orientation
     vxy_ang = normalize_angle(np.arctan2(vy, vx), lower=0, upper=np.pi, dtype=np.float32)