diff --git a/MRISnapshot/create_report.py b/MRISnapshot/create_report.py
index 90e32f7..310fb98 100644
--- a/MRISnapshot/create_report.py
+++ b/MRISnapshot/create_report.py
@@ -52,12 +52,14 @@ def parse_config(df_conf, list_col_names):
     cols_default = ['id_col', 'ulay_col', 'mask_col', 'olay_col', 'olay_col2', 'sel_vals_olay', 
                     'sel_vals_olay2', 'view_plane', 'num_slice', 'step_size_slice',
                     'min_vox', 'crop_to_mask', 'crop_to_olay', 'padding_ratio', 'bin_olay', 
+                    'segment_olay', 'num_classes_olay',
                     'is_edge', 'alpha_olay', 'perc_high', 'perc_low', 
                     'is_out_single', 'is_out_noqc', 'img_width',
                     'label_checkbox1', 'label_checkbox2', 'label_editbox']
     vals_default = ['', '', '', '', '', '', 
                     '', 'A', '4', '',
-                    '1', '0', '0', '', '1', 
+                    '1', '0', '0', '', '1',
+                    '0', '0',
                     '1', '1', '99', '1', 
                     '1', '0', '300',
                     'PASS', 'FAIL', 'Notes']
@@ -104,7 +106,8 @@ def parse_config(df_conf, list_col_names):
     params.sel_vals_olay2 = [int(n) for n in params.sel_vals_olay2.split('+') if n != '']
 
     ### Convert numeric args from str to int or float
-    for tmp_arg in ['num_slice', 'step_size_slice', 'min_vox', 'crop_to_mask', 'crop_to_olay', 'bin_olay', 
+    for tmp_arg in ['num_slice', 'step_size_slice', 'min_vox', 'crop_to_mask', 'crop_to_olay', 'bin_olay',
+                    'segment_olay', 'num_classes_olay',
                     'is_edge', 'is_out_single', 'is_out_noqc', 'img_width']:
         if params[tmp_arg] != '':
             params[tmp_arg] = int(params[tmp_arg])
@@ -118,6 +121,20 @@ def parse_config(df_conf, list_col_names):
         params.alpha_olay = 1
         logger.warning('    Parameter alpha_olay is set to 1, because is_edge was set to 1')
 
+    ### Correct inconsistent parameters
+    if (params.segment_olay == 1):
+        if params.num_classes_olay > 10:
+            params.num_classes_olay = 10
+            logger.warning('    Parameter num_classes_olay is set to 10 (max value)')
+        if params.num_classes_olay < 2:
+            params.num_classes_olay = 2
+            logger.warning('    Parameter num_classes_olay is set to 2 (min value)')
+        if params.is_edge == 0:
+            params.is_edge = 1
+            logger.warning('    Parameter is_edge is set to 1, because segment_olay was set to 1')
+        params.alpha_olay = 1
+        logger.warning('    Parameter alpha_olay is set to 1, because is_edge was set to 1')
+
     return params
 
 def create_dir(dir_name):
@@ -199,6 +216,14 @@ def read_and_check_images(df_images, params, sub_index, orient = 'LPS'):
             else:
                 try:
                     nii = nib.load(fname)
+                except:
+                    nii = None
+                    qc_ok_flag = 0
+                    qc_msg = 'Could not read ' + col_name + ' image ' + fname
+                    logger.warning('   ' + qc_msg + ', subject discarded!')
+                    return qc_ok_flag, qc_msg, nii_out, fnames_out
+
+                try:
                     orig_ornt = nib.io_orientation(nii.affine)
                     targ_ornt = axcodes2ornt(orient)
                     if np.all(orig_ornt == targ_ornt) == False:
@@ -207,17 +232,22 @@ def read_and_check_images(df_images, params, sub_index, orient = 'LPS'):
                     if ref_affine == []:
                         ref_affine = nii.affine
                     else:
-                        if (ref_affine - nii.affine).sum() != 0:
+
+                        ## We compare the affine matrix for underlay and overlay images
+                        ##   - we use 3 digit precision
+                        ##     otherwise it may be too sensitive to numeric differences in headers
+                        if (np.abs(ref_affine - nii.affine).sum()) > 0.001:
                             qc_ok_flag = 0
-                            qc_msg = 'Inconsistent image vs. overlay'
+                            qc_msg = 'Inconsistent image vs. overlay: affine matrices in headers not consistent'
                             logger.warning('   ' + qc_msg + ', subject discarded!')
                             return qc_ok_flag, qc_msg, nii_out, fnames_out
                 except:
                     nii = None
                     qc_ok_flag = 0
-                    qc_msg = 'Could not read ' + col_name + ' image ' + fname
+                    qc_msg = 'Could not reorient ' + col_name + ' image ' + fname
                     logger.warning('   ' + qc_msg + ', subject discarded!')
                     return qc_ok_flag, qc_msg, nii_out, fnames_out
+                    
                         
         nii_out.append(nii)
         fnames_out.append(fname)
@@ -322,6 +352,37 @@ def scale_img_contrast(nii_img, nii_mask, perc_low, perc_high):
 ']')
     return nii_out
 
+def digitize_olay(nii, num_classes, perc_low = 5, perc_high = 95):
+    '''Simple quick segmentation of overlay image by quantizing intensities.
+    The purpose is to make sure that edge detection on overlay will be smooth. 
+    Bins for digitization will be estimated from robust intensity ranges
+    (using percentile values to exclude outliers)
+    
+    :param nii: Input nifti image
+    :param num_classes: Number of classes (digits) in output image
+    
+    :return nii_out: Output nifti image
+    '''
+    ## Check input img
+    if nii == None:
+        return nii
+    
+    ## Get img data
+    tmp_img = nii.get_fdata()
+    
+    ## Find min and max intensities
+    pmin, pmax = np.percentile(tmp_img.flatten(), [perc_low, perc_high])
+    
+    ## Digitize image
+    tmp_img = np.digitize(tmp_img, np.arange(pmin, pmax, (pmax-pmin)/num_classes))
+    
+    ## Create out nifti
+    nii_out = nib.Nifti1Image(tmp_img, nii.affine, nii.header)
+
+    ## Return out nifti
+    logger.info('      Overlay image segmented to ' + str(num_classes) + ' classes')
+    return nii_out
+
 def extract_snapshot(img_ulay, img_olay, img_olay2, params, curr_view, curr_slice, slice_index, 
                      sub_id, dir_snapshots_full, list_sel_slices):
     ''' Extracts an image snapshot based on input parameters, and writes it to 
@@ -345,7 +406,8 @@ def extract_snapshot(img_ulay, img_olay, img_olay2, params, curr_view, curr_slic
     img2d_ulay = img_ulay[:,:,curr_slice].astype(float)
     
     # Scale underlay image between 0 and 1
-    img2d_ulay = (img2d_ulay - img2d_ulay.min()) / (img2d_ulay.max() - img2d_ulay.min())
+    if img2d_ulay.max() - img2d_ulay.min() > 0:
+        img2d_ulay = (img2d_ulay - img2d_ulay.min()) / (img2d_ulay.max() - img2d_ulay.min())
 
     ## Resize underlay slice
     #img2d_ulay = zoom(img2d_ulay, (scX,scY), order=1)
@@ -633,6 +695,11 @@ def create_snapshots(params, df_images, dir_snapshots_full, out_dir):
                     ## Scale ulay image intensities
                     nii_ulay = scale_img_contrast(nii_ulay, nii_mask, params.perc_low, params.perc_high)
 
+                    ## Digitize olay image intensities
+                    if params.segment_olay == 1:
+                        nii_olay = digitize_olay(nii_olay, params.num_classes_olay)
+                        nii_olay2 = digitize_olay(nii_olay2, params.num_classes_olay)
+
                     ### Create snapshots for each orientation
                     for view_index, curr_view in enumerate(params.view_plane):
 
diff --git a/MRISnapshot/prep_data.py b/MRISnapshot/prep_data.py
index a2a2284..a2ca74f 100644
--- a/MRISnapshot/prep_data.py
+++ b/MRISnapshot/prep_data.py
@@ -115,8 +115,10 @@ def prep_dataset(params):
                         'crop_to_mask' : 0,
                         'crop_to_olay' : 0,
                         'padding_ratio' : 0,
-                        'bin_olay' : 0,  
-                        'is_edge' : 1, 
+                        'bin_olay' : 0,
+                        'segment_olay' : 0,
+                        'num_classes_olay' : 0,
+                        'is_edge' : 1,
                         'alpha_olay' : 1, 
                         'perc_high' : 100, 'perc_low' : 0, 
                         'is_out_single' : 0, 'is_out_noqc' : 0,