diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..b3ada6a
Binary files /dev/null and b/.DS_Store differ
diff --git a/analysis/analysis.py b/analysis/analysis.py
index bba0edd..30b495a 100755
--- a/analysis/analysis.py
+++ b/analysis/analysis.py
@@ -23,12 +23,12 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 # matplotlib controls
 plt.rcParams['svg.fonttype'] = 'none'  # to embed fonts in output ('path' is to convert as text as paths)
@@ -39,7 +39,7 @@
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # dimensions
     params['dimensions']={}
@@ -66,20 +66,16 @@ def hist_dimensions(cells, outputdir='.', bins=['auto','auto','auto','auto'], un
 
     # initialization
     if mode == 'um':
-        titles = [u'length (\u03BCm)', u'width (\u03BCm)', u'area (\u03BCm\u00B2)', u'volume (\u03BCm\u00B3)']
+        titles = ['length (\u03BCm)', 'width (\u03BCm)', 'area (\u03BCm\u00B2)', 'volume (\u03BCm\u00B3)']
         attrs = ['height_um','width_um','area_um2', 'volume_um3']
     else:
-        titles = [u'length (px)', u'width (px)', u'area (px\u00B2)', u'volume (px\u00B3)']
+        titles = ['length (px)', 'width (px)', 'area (px\u00B2)', 'volume (px\u00B3)']
         attrs = ['height','width','area', 'volume']
 
     nattrs = len(attrs)
-    if bins is None:
-        bins = ['auto' for i in range(nattrs)]
-    elif len(bins) != nattrs:
+    if len(bins) != nattrs:
         raise ValueError("bins has the wrong dimensions!")
-    if units_dx is None:
-        units_dx = [None for i in range(nattrs)]
-    elif len(units_dx) != nattrs:
+    if len(units_dx) != nattrs:
         raise ValueError("units_dx has the wrong dimensions!")
     ncells = len(cells)
 
@@ -90,7 +86,7 @@ def hist_dimensions(cells, outputdir='.', bins=['auto','auto','auto','auto'], un
 
     # make lists
     data = [ [] for i in range(nattrs)]
-    keys = cells.keys()
+    keys = list(cells.keys())
     for n in range(ncells):
         key = keys[n]
         cell = cells[key]
@@ -118,13 +114,13 @@ def hist_dimensions(cells, outputdir='.', bins=['auto','auto','auto','auto'], un
 
     for i in range(nattrs):
         attr = attrs[i]
-        print "attr = {}".format(attr)
+        print("attr = {}".format(attr))
         ax = axes[i]
 
         # compute histogram
         hist,edges = np.histogram(data[i], bins=bins[i], density=False)
         nbins = len(edges)-1
-        print "nbins = {:d}".format(nbins)
+        print("nbins = {:d}".format(nbins))
 
         # plot histogram
         color = 'grey'
@@ -165,7 +161,7 @@ def hist_dimensions(cells, outputdir='.', bins=['auto','auto','auto','auto'], un
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=0.):
@@ -203,7 +199,7 @@ def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=
 
     # make lists
     data = [ [] for i in range(nattrs)]
-    keys = cells.keys()
+    keys = list(cells.keys())
     for n in range(ncells):
         key = keys[n]
         cell = cells[key]
@@ -230,7 +226,7 @@ def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=
         axes.append(ax)
 
     for i in range(nattrs):
-        print "attr number = {:d}".format(i)
+        print("attr number = {:d}".format(i))
         ax = axes[i]
 
         # remove outliers
@@ -238,7 +234,7 @@ def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=
         ntot = len(data_sorted)
         n0 = int(0.5*alpha*float(ntot))
         n1 = ntot - n0
-        print "{:2s}n0 = {:d}    n1 = {:d}".format("", n0, n1)
+        print("{:2s}n0 = {:d}    n1 = {:d}".format("", n0, n1))
         data_trunc = data_sorted[n0:n1]
 
         mu = np.nanmean(data_trunc)
@@ -250,7 +246,7 @@ def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=
         # compute histogram
         hist,edges = np.histogram(data_trunc, bins=bins[i], density=False)
         nbins = len(edges)-1
-        print "{:2s}nbins = {:d}".format("",nbins)
+        print("{:2s}nbins = {:d}".format("",nbins))
 
         # plot histogram
         color = 'grey'
@@ -281,7 +277,7 @@ def hist_dimensions_other(cells, outputdir='.', bins=None, units_dx=None, alpha=
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, mode='total',qcut=0, bg_val=None, backgrounds=None):
@@ -291,11 +287,11 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
 
     # initialization
     if mode == 'concentration_fl':
-        print "concentration_fl mode"
+        print("concentration_fl mode")
         fl_dtype = np.float_
         fmt = "$\\mu = {mu:,.2f}$\n$\\sigma = {sig:,.2f}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,.2f}$"
     elif mode == 'total_fl':
-        print "total_fl mode"
+        print("total_fl mode")
         fl_dtype = np.uint16
         fmt = "$\\mu = {mu:,d}$\n$\\sigma = {sig:,d}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,d}$"
     else:
@@ -304,7 +300,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
     data_bg = []
 
     # filling up the fluorescence
-    cellref = cells.values()[0]
+    cellref = list(cells.values())[0]
     nchannel = len(cellref['fluorescence']['total'])
 
     if bins is None:
@@ -322,7 +318,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
         BG = []
 
         # make lists
-        keys = cells.keys()
+        keys = list(cells.keys())
         for n in range(ncells):
             key = keys[n]
             cell = cells[key]
@@ -360,7 +356,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
     else:
         bgcolor='g'
         bgs = [np.float_(backgrounds[i]) for i in range(nchannel)]
-        #print bgs
+        print(bgs)
 
     # make figure
     fig = plt.figure(num=None, facecolor='w', figsize=(nchannel*4,3))
@@ -373,7 +369,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
         axes.append(ax)
 
     for i in range(nchannel):
-        print "channel = {}".format(i)
+        print("channel = {}".format(i))
         ax = axes[i]
 
         # compute histogram
@@ -385,7 +381,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
         d = d[n0:n1+1]
         hist,edges = np.histogram(d, bins=bins[i], density=False)
         nbins = len(edges)-1
-        print "nbins = {:d}".format(nbins)
+        print("nbins = {:d}".format(nbins))
 
         # plot histogram
         color = 'grey'
@@ -434,7 +430,7 @@ def hist_channels(cells, outputdir='.', bins=None, units_dx=None, titles=None, m
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['I1','I2','QUEEN'], channels=[0,1], colors=['darkblue', 'darkgreen', 'darkblue'], units_dx=[None, None, None], mode='total_fl', bgcolor='r',qcut=0, backgrounds=None):
@@ -455,12 +451,12 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     QUEEN=[]
 
     if mode == 'concentration_fl':
-        print "concentration_fl mode"
+        print("concentration_fl mode")
         fl_dtype = np.float_
         fl_fmt = "$\\mu = {mu:,.2f}$\n$\\sigma = {sig:,.2f}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,.2f}$"
 
     elif mode == 'total_fl':
-        print "total_fl mode"
+        print("total_fl mode")
         fl_dtype = np.uint16
         fl_fmt = "$\\mu = {mu:,d}$\n$\\sigma = {sig:,d}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,d}$"
     else:
@@ -472,7 +468,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         bgcolor='g'
 
     # make lists
-    keys = cells.keys()
+    keys = list(cells.keys())
     for n in range(ncells):
         key = keys[n]
         cell = cells[key]
@@ -498,12 +494,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
                 bg_y  = backgrounds[c2]
         elif mode == 'total_fl':
             pass
-        try:
-            z = (float(x)-float(bg_x))/(float(y)-float(bg_y))
-        except ZeroDivisionError:
-            print("x = {:.6f}    bg_x = {:.6f}    y = {:.6f}    bg_y = {:.6f}".format(x,bg_x,y,bg_y))
-            print("Skipping cell {:s}".format(cell['id']))
-            continue
+        z = (float(x)-float(bg_x))/(float(y)-float(bg_y))
         I1.append(x)
         BG1.append(bg_x)
         I2.append(y)
@@ -516,8 +507,8 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     I2 = np.array(I2, dtype=fl_dtype)
     BG2 = np.array(BG2, dtype=fl_dtype)
     QUEEN = np.array(QUEEN, dtype=np.float_)
-    #print np.unique(BG1)
-    #print np.unique(BG2)
+    print(np.unique(BG1))
+    print(np.unique(BG2))
 
     # make plot
     N = len(I1)
@@ -529,7 +520,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     mu_bg1 = np.median(BG1)
     mu_bg2 = np.median(BG2)
     data_bg = [BG1, BG2]
-    print "bg1 = {:.2f}    bg2 = {:.2f}".format(mu_bg1,mu_bg2)
+    print("bg1 = {:.2f}    bg2 = {:.2f}".format(mu_bg1,mu_bg2))
     fmts = [fl_fmt,fl_fmt,"$\\mu = {mu:.2f}$\n$\\sigma = {sig:.2f}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:.2f}$"]
 
     fig = plt.figure(num=None, facecolor='w', figsize=(3*4,3))
@@ -553,7 +544,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         d = d[n0:n1+1]
         hist,edges = np.histogram(d, bins=bins[i], density=False)
         nbins = len(edges)-1
-        print "nbins = {:d}".format(nbins)
+        print("nbins = {:d}".format(nbins))
 
         # plot histogram
         color = colors[i]
@@ -606,7 +597,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 #################### main ####################
@@ -628,7 +619,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
 
     cells = load_json2dict(cellfile)
     ncells = len(cells)
-    print "ncells = {:d}".format(ncells)
+    print("ncells = {:d}".format(ncells))
 
     # output directory
     outputdir = namespace.outputdir
@@ -639,7 +630,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     outputdir = os.path.join(outputdir,'analysis')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
@@ -660,7 +651,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         mydir = os.path.join(outputdir,'dimensions')
         if not os.path.isdir(mydir):
             os.makedirs(mydir)
-        print "{:<20s}{:<s}".format("outputdir", mydir)
+        print("{:<20s}{:<s}".format("outputdir", mydir))
         hist_dimensions(cells, outputdir=mydir, **allparams['dimensions'])
 
     # dimensions other
@@ -668,7 +659,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         mydir = os.path.join(outputdir,'dimensions')
         if not os.path.isdir(mydir):
             os.makedirs(mydir)
-        print "{:<20s}{:<s}".format("outputdir", mydir)
+        print("{:<20s}{:<s}".format("outputdir", mydir))
         hist_dimensions_other(cells, outputdir=mydir, **allparams['dimensions_other'])
 
     # make queen analysis
@@ -676,7 +667,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         mydir = os.path.join(outputdir,'fluorescence')
         if not os.path.isdir(mydir):
             os.makedirs(mydir)
-        print "{:<20s}{:<s}".format("outputdir", mydir)
+        print("{:<20s}{:<s}".format("outputdir", mydir))
         hist_channels(cells, outputdir=mydir, **allparams['fluorescence'])
 
     # make queen analysis
@@ -684,6 +675,6 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
         mydir = os.path.join(outputdir,'queen')
         if not os.path.isdir(mydir):
             os.makedirs(mydir)
-        print "{:<20s}{:<s}".format("outputdir", mydir)
+        print("{:<20s}{:<s}".format("outputdir", mydir))
         make_plot_queen(cells, outputdir=mydir, **allparams['queen'])
 
diff --git a/analysis/analysis_overlays.py b/analysis/analysis_overlays.py
index 1e21850..35d0368 100755
--- a/analysis/analysis_overlays.py
+++ b/analysis/analysis_overlays.py
@@ -13,7 +13,6 @@
 import matplotlib.pyplot as plt
 import matplotlib.gridspec as mgs
 import matplotlib.ticker
-import matplotlib.patches as mpatches
 from matplotlib.colors import cnames as colornames
 
 # custom
@@ -25,12 +24,12 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 # matplotlib controls
 plt.rcParams['svg.fonttype'] = 'none'  # to embed fonts in output ('path' is to convert as text as paths)
@@ -41,7 +40,7 @@
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['analysis_overlays']={}
@@ -63,17 +62,17 @@ def default_parameters():
     mydict['mode']='total_fl'  # alternative value is \'concentration\'
     return params
 
-def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, units_dx=None, title=None, mode='um',qcut=0, ncol=None):
+def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, units_dx=None, title=None, mode='um',qcut=0):
     """
     Make an histogram of the signal obtained per cell.
     """
 
     # initialization
     if mode == 'um':
-        titles = [u'length (\u03BCm)', u'width (\u03BCm)', u'area (\u03BCm\u00B2)', u'volume (\u03BCm\u00B3)']
+        titles = ['length (\u03BCm)', 'width (\u03BCm)', 'area (\u03BCm\u00B2)', 'volume (\u03BCm\u00B3)']
         attrs = ['height_um','width_um','area_um2', 'volume_um3']
     else:
-        titles = [u'length (px)', u'width (px)', u'area (px\u00B2)', u'volume (px\u00B3)']
+        titles = ['length (px)', 'width (px)', 'area (px\u00B2)', 'volume (px\u00B3)']
         attrs = ['height','width','area', 'volume']
 
     nattrs = len(attrs)
@@ -83,11 +82,11 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
     ndata = len(celldicts)
     if ndata == 0:
         raise ValueError("Empty input data")
-    print "ndata = {:d}".format(ndata)
+    print("ndata = {:d}".format(ndata))
 
     # colors
     if colors is None:
-        colors = colornames.keys()[:ndata]
+        colors = list(colornames.keys())[:ndata]
 
     # bins
     if bins is None:
@@ -114,7 +113,7 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
 
         # make lists
         dimensions = [ [] for i in range(nattrs)]
-        keys = cells.keys()
+        keys = list(cells.keys())
         for n in range(ncells):
             key = keys[n]
             cell = cells[key]
@@ -146,11 +145,11 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
 
     for j in range(nattrs):
         attr = attrs[j]
-        print "attr {:d} / {:d}".format(j,nattrs-1)
+        print("attr {:d} / {:d}".format(j,nattrs-1))
         ax = axes[j]
 
         for i in range(ndata):
-            print "{:<2s}data {:d} / {:d}".format("", i, ndata-1)
+            print("{:<2s}data {:d} / {:d}".format("", i, ndata-1))
             # compute histogram
             d = data[i][j]
             N = len(d)
@@ -158,14 +157,18 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
             n0 = int(qcut*float(N))
             n1 = min(int((1.-qcut)*float(N)),N-1)
             d = d[n0:n1+1]
-            print "{:<4s}qcut = {:.1f} %".format("",qcut*100)
+            print("{:<4s}qcut = {:.1f} %".format("",qcut*100))
             hist,edges = np.histogram(d, bins=bins[j][i], density=True)
-            print "{:<4s}nbins = {:,d}".format("",len(edges)-1)
+            print("{:<4s}nbins = {:,d}".format("",len(edges)-1))
 
             # plot histogram
             color = colors[i]
             #ax.bar(edges[:-1], hist, np.diff(edges), color='none', edgecolor=color, lw=0.5, label=labels[i])
-            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5)
+            if j == 0:
+                label = labels[i]
+            else:
+                label = None
+            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5, label=label)
         # end loop
     # end loop
 
@@ -176,8 +179,8 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
         #ax.annotate(fmts[i].format(mu=mus[i],sig=sigs[i], N=len(data[i]), med=meds[i]), xy=(0.70,0.98), xycoords='axes fraction', ha='left', va='top')
 
         # adjust the axis
-#        if (j == 0):
-#            ax.legend(loc='best',fontsize="medium",frameon=False)
+        if (j == 0):
+            ax.legend(loc='best',fontsize="medium",frameon=False)
         #ax.set_ylabel("pdf",fontsize="medium",labelpad=10)
         ax.tick_params(length=4)
         ax.tick_params(axis='both', labelsize='medium', labelleft=False, left=False)
@@ -192,25 +195,7 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
         ax.spines['right'].set_visible(False)
         ax.spines['top'].set_visible(False)
 
-    # legend
-    patches=[]
-    for n in range(ndata):
-        label=labels[n]
-        color=colors[n]
-        if label is None:
-            label = "{:d}".format(n)
-        patch = mpatches.Patch(color=color, label=label)
-        patches.append(patch)
-    #axes[2].add_artist(plt.legend(handles=patches, loc='upper right'))
-    if (ncol is None):
-        ncol = ndata
-    #plt.figlegend(handles=patches, fontsize='small', mode='expand', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-    plt.figlegend(handles=patches, fontsize='small', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-
-    nrow=int(float(ndata)/ncol + 0.5)
-    a = 0.05
-    rect=[0., nrow*a, 1., 1.]
-    gs.tight_layout(fig, rect=rect)
+    gs.tight_layout(fig)
     if mode == 'um':
         filename = 'analysis_overlay_dimensions_um'
     else:
@@ -220,21 +205,21 @@ def hist_dimensions(celldicts, labels, outputdir='.', bins=None, colors=None, un
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
-def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units_dx=None, titles=None, mode='total_fl',qcut=0, xminxmax=None, backgrounds=None, ncol=None):
+def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units_dx=None, titles=None, mode='total_fl',qcut=0, xminxmax=None, backgrounds=None):
     """
     Make an histogram of the signal obtained per cell.
     """
 
     # initialization
     if mode == 'concentration_fl':
-        print "concentration_fl mode"
+        print("concentration_fl mode")
         fl_dtype = np.float_
         fmt = "$\\mu = {mu:,.2f}$\n$\\sigma = {sig:,.2f}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,.2f}$"
     elif mode == 'total_fl':
-        print "total_fl mode"
+        print("total_fl mode")
         fl_dtype = np.uint16
         fmt = "$\\mu = {mu:,d}$\n$\\sigma = {sig:,d}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,d}$"
     else:
@@ -246,15 +231,15 @@ def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units
         raise ValueError("Empty input data")
     if ndata != len(labels):
         raise ValueError("Labels should have the same dimension as cell dicts")
-    print "ndata = {:d}".format(ndata)
+    print("ndata = {:d}".format(ndata))
 
     # channels
-    nchannel = len(celldicts[0].values()[0]['fluorescence']['total'])
-    print "nchannel = {:d}".format(nchannel)
+    nchannel = len(list(celldicts[0].values())[0]['fluorescence']['total'])
+    print("nchannel = {:d}".format(nchannel))
 
     # colors
     if colors is None:
-        colors = colornames.keys()[:ndata]
+        colors = list(colornames.keys())[:ndata]
 
     # bins
     if bins is None:
@@ -288,7 +273,7 @@ def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units
             BG = []
 
             # make lists
-            keys = cells.keys()
+            keys = list(cells.keys())
             for n in range(ncells):
                 key = keys[n]
                 cell = cells[key]
@@ -339,26 +324,31 @@ def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units
         axes.append(ax)
 
     for c in range(nchannel):
-        print "channel = {:d} / {:d}".format(c,nchannel-1)
+        print("channel = {:d} / {:d}".format(c,nchannel-1))
         ax=axes[c]
         for i in range(ndata):
-            print "{:<2s}data {:d} / {:d}".format("", i, ndata-1)
+            print("{:<2s}data {:d} / {:d}".format("", i, ndata-1))
             # compute histogram
             d = data[c][i]
             N = len(d)
             d = np.sort(d)
             n0 = int(0.5*qcut*float(N))
-            print "{:<4s}qcut = {:.1f} %".format("",qcut*100)
+            print("{:<4s}qcut = {:.1f} %".format("",qcut*100))
             n1 = N - n0
 #            print n0, n1
             d = d[n0:n1]
             hist,edges = np.histogram(d, bins=bins[c][i], density=True)
-            print "{:<4s}nbins = {:,d}".format("",len(edges)-1)
+            print("{:<4s}nbins = {:,d}".format("",len(edges)-1))
 
             # plot histogram
             color = colors[i]
             #ax.bar(edges[:-1], hist, np.diff(edges), color='none', edgecolor=color, lw=0.5, label=labels[i])
-            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5)
+            if (c ==0):
+                label=labels[i]
+            else:
+                label=None
+
+            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5, label=label)
             # end loop on data
 
         # plot background
@@ -387,25 +377,7 @@ def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units
         ax.spines['top'].set_visible(False)
     # end loop on channels
 
-    # legend
-    patches=[]
-    for n in range(ndata):
-        label=labels[n]
-        color=colors[n]
-        if label is None:
-            label = "{:d}".format(n)
-        patch = mpatches.Patch(color=color, label=label)
-        patches.append(patch)
-    #axes[2].add_artist(plt.legend(handles=patches, loc='upper right'))
-    if (ncol is None):
-        ncol = ndata
-    #plt.figlegend(handles=patches, fontsize='small', mode='expand', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-    plt.figlegend(handles=patches, fontsize='small', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-
-    nrow=int(float(ndata)/ncol + 0.5)
-    a = 0.05
-    rect=[0., nrow*a, 1., 1.]
-    gs.tight_layout(fig, rect=rect)
+    gs.tight_layout(fig)
     if mode == 'concentration_fl':
         filename = 'analysis_overlay_concentration_fl'
     elif mode == 'total_fl':
@@ -415,21 +387,21 @@ def hist_channel(celldicts, labels, outputdir='.', bins=None, colors=None, units
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
-def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colors=None, units_dx=None, titles=None, mode='total_fl',qcut=0, xminxmax=None, backgrounds=None, ncol=None):
+def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colors=None, units_dx=None, titles=None, mode='total_fl',qcut=0, xminxmax=None, backgrounds=None):
     """
     Make an histogram of the signal obtained per cell for the fluorescence channels used to measure the Queen ratio. Also plot the distribution of the Queen ratio.
     """
 
     # initialization
     if mode == 'concentration_fl':
-        print "concentration_fl mode"
+        print("concentration_fl mode")
         fl_dtype = np.float_
         fmt = "$\\mu = {mu:,.2f}$\n$\\sigma = {sig:,.2f}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,.2f}$"
     elif mode == 'total_fl':
-        print "total_fl mode"
+        print("total_fl mode")
         fl_dtype = np.uint16
         fmt = "$\\mu = {mu:,d}$\n$\\sigma = {sig:,d}$\n$N = {N:,d}$\n$\\mathrm{{med}} = {med:,d}$"
     else:
@@ -441,18 +413,18 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
         raise ValueError("Empty input data")
     if ndata != len(labels):
         raise ValueError("Labels should have the same dimension as cell dicts")
-    print "ndata = {:d}".format(ndata)
+    print("ndata = {:d}".format(ndata))
 
     # channels
     nchannel = 2
     c1 = channels[0]
     c2 = channels[1]
-    print "nchannel = {:d}, c1 = {:d}, c2 = {:d}".format(nchannel, c1, c2)
+    print("nchannel = {:d}, c1 = {:d}, c2 = {:d}".format(nchannel, c1, c2))
     nplot = nchannel + 1
 
     # colors
     if colors is None:
-        colors = colornames.keys()[:ndata]
+        colors = list(colornames.keys())[:ndata]
 
     # bins
     if bins is None:
@@ -491,7 +463,7 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
         QUEEN = []
 
         # make lists
-        keys = cells.keys()
+        keys = list(cells.keys())
         for n in range(ncells):
             key = keys[n]
             cell = cells[key]
@@ -518,12 +490,7 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
             elif mode == 'total_fl':
                 sys.exit("Not programmed yet!")
                 pass
-            try:
-                z = ((float(x1) -  float(bg_x1)) / (float(x2) - float(bg_x2)) )
-            except ZeroDivisionError:
-                print("x = {:.6f}    bg_x = {:.6f}    y = {:.6f}    bg_y = {:.6f}".format(x1,bg_x1,x2,bg_x2))
-                print("Skipping cell {:s}".format(cell['id']))
-                continue
+            z = ((float(x1) -  float(bg_x1)) / (float(x2) - float(bg_x2)) )
             FL1.append(x1)
             FL2.append(x2)
             BG1.append(bg_x1)
@@ -557,45 +524,52 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
     # plot fluorescence
     data=[data_I1, data_I2]
     for c in [c1,c2]:
-        print "channel = {:d} / {:d}".format(c,1)
+        print("channel = {:d} / {:d}".format(c,1))
         ax=axes[c]
         for i in range(ndata):
-            print "{:<2s}data {:d} / {:d}".format("", i, ndata-1)
+            print("{:<2s}data {:d} / {:d}".format("", i, ndata-1))
             # compute histogram
             d = data[c][i]
             N = len(d)
             d = np.sort(d)
-            n0 = int(qcut*float(N))
-            n1 = min(int((1.-qcut)*float(N)),N-1)
-            d = d[n0:n1+1]
+            n0 = int(0.5*qcut*float(N))
+            print("{:<4s}qcut = {:.1f} %".format("",qcut*100))
+            n1 = N - n0
+#            print n0, n1
+            d = d[n0:n1]
             hist,edges = np.histogram(d, bins=bins[c][i], density=True)
-            nbins = len(edges)-1
-            print "{:<4s}nbins = {:,d}".format("",nbins)
+            print("{:<4s}nbins = {:,d}".format("",len(edges)-1))
 
             # plot histogram
             color = colors[i]
             #ax.bar(edges[:-1], hist, np.diff(edges), color='none', edgecolor=color, lw=0.5, label=labels[i])
-            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5)
+            if (c ==0):
+                label=labels[i]
+            else:
+                label=None
+
+            ax.plot(0.5*(edges[:-1]+edges[1:]), hist, '-', color=color, lw=0.5, label=label)
 
         # plot background
         ax.axvline(x=bgs[c], color=bgcolor, lw=0.5, ls='--')
-        print "background = {:.2f}".format(bgs[c])
+        print("background = {:.2f}".format(bgs[c]))
         # end loop on data
 
     # plot Queen ratio
     ax=axes[2]
     for i in range(ndata):
-        print "{:<2s}data {:d} / {:d}".format("", i, ndata-1)
+        print("{:<2s}data {:d} / {:d}".format("", i, ndata-1))
         # compute histogram
         d = data_queen[i]
         N = len(d)
         d = np.sort(d)
-        n0 = int(qcut*float(N))
-        n1 = min(int((1.-qcut)*float(N)),N-1)
-        d = d[n0:n1+1]
+        n0 = int(0.5*qcut*float(N))
+        print("{:<4s}qcut = {:.1f} %".format("",qcut*100))
+        n1 = N - n0
+#            print n0, n1
+        d = d[n0:n1]
         hist,edges = np.histogram(d, bins=bins[2][i], density=True)
-        nbins = len(edges)-1
-        print "{:<4s}nbins = {:,d}".format("",nbins)
+        print("{:<4s}nbins = {:,d}".format("",len(edges)-1))
 
         # plot histogram
         color = colors[i]
@@ -633,25 +607,7 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
         ax.spines['top'].set_visible(False)
     # end loop on channels
 
-    # legend
-    patches=[]
-    for n in range(ndata):
-        label=labels[n]
-        color=colors[n]
-        if label is None:
-            label = "{:d}".format(n)
-        patch = mpatches.Patch(color=color, label=label)
-        patches.append(patch)
-    #axes[2].add_artist(plt.legend(handles=patches, loc='upper right'))
-    if (ncol is None):
-        ncol = ndata
-    #plt.figlegend(handles=patches, fontsize='small', mode='expand', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-    plt.figlegend(handles=patches, fontsize='small', ncol=ncol,borderaxespad=0, borderpad=0, loc='lower center', frameon=False)
-
-    nrow=int(float(ndata)/ncol + 0.5)
-    a = 0.05
-    rect=[0., nrow*a, 1., 1.]
-    gs.tight_layout(fig, rect=rect)
+    gs.tight_layout(fig)
     if mode == 'concentration_fl':
         filename = 'analysis_overlay_queen_concentration_fl'
     elif mode == 'total_fl':
@@ -661,7 +617,7 @@ def hist_queen(celldicts, labels, outputdir='.', channels=[0,1], bins=None, colo
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto'], channels=[0,1], colors=['darkblue', 'darkgreen', 'darkblue'], units_dx=None, qcut=0):
@@ -679,7 +635,7 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
     ndata = len(celldicts)
     if ndata == 0:
         raise ValueError("Empty input data")
-    print "ndata = {:d}".format(ndata)
+    print("ndata = {:d}".format(ndata))
 
     for i in range(ndata):
         cells = celldicts[i]
@@ -691,7 +647,7 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
         QUEEN=[]
 
         # make lists
-        keys = cells.keys()
+        keys = list(cells.keys())
         for n in range(ncells):
             key = keys[n]
             cell = cells[key]
@@ -732,7 +688,7 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
         n1 = min(int((1.-qcut)*float(N)),N-1)
         d = d[n0:n1+1]
         hist,edges = np.histogram(d, bins=bins[i], density=True)
-        print "nbins = {:,d}".format(len(edges)-1)
+        print("nbins = {:,d}".format(len(edges)-1))
 
         # plot histogram
         color = colors[i]
@@ -773,7 +729,7 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 #################### main ####################
@@ -800,10 +756,10 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
     celldicts=[]
     for n in range(nfiles):
         cellfile = cellfiles[n]
-        print cellfile
+        print(cellfile)
         cells = load_json2dict(cellfile)
         ncells = len(cells)
-        print "ncells = {:d}".format(ncells)
+        print("ncells = {:d}".format(ncells))
         celldicts.append(cells)
 
     # labels
@@ -824,7 +780,7 @@ def hist_queen_old(celldicts, labels, outputdir='.', bins=['auto','auto','auto']
     outputdir = os.path.join(rootdir,'overlays')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
diff --git a/analysis/analysis_queen.py b/analysis/analysis_queen.py
index 92ffd18..6f3ac1b 100755
--- a/analysis/analysis_queen.py
+++ b/analysis/analysis_queen.py
@@ -23,12 +23,12 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 # matplotlib controls
 plt.rcParams['svg.fonttype'] = 'none'  # to embed fonts in output ('path' is to convert as text as paths)
@@ -39,7 +39,7 @@
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['analysis_queen']={}
@@ -68,7 +68,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     QUEEN=[]
 
     # make lists
-    keys = cells.keys()
+    keys = list(cells.keys())
     for n in range(ncells):
         key = keys[n]
         cell = cells[key]
@@ -145,7 +145,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     for ext in exts:
         fileout = os.path.join(outputdir,filename+ext)
         fig.savefig(fileout, bbox_inches='tight', pad_inches=0)
-        print "Fileout: {:<s}".format(fileout)
+        print("Fileout: {:<s}".format(fileout))
     return
 
 #################### main ####################
@@ -167,7 +167,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
 
     cells = load_json2dict(cellfile)
     ncells = len(cells)
-    print "ncells = {:d}".format(ncells)
+    print("ncells = {:d}".format(ncells))
 
     # output directory
     outputdir = namespace.outputdir
@@ -179,7 +179,7 @@ def make_plot_queen(cells, outputdir='.', bins=['auto','auto','auto'], titles=['
     outputdir = os.path.join(rootdir,'queen')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
diff --git a/dockerfile/.DS_Store b/dockerfile/.DS_Store
new file mode 100644
index 0000000..24da4e8
Binary files /dev/null and b/dockerfile/.DS_Store differ
diff --git a/image_processing/__pycache__/acquisitionlib.cpython-36.pyc b/image_processing/__pycache__/acquisitionlib.cpython-36.pyc
new file mode 100644
index 0000000..076a16f
Binary files /dev/null and b/image_processing/__pycache__/acquisitionlib.cpython-36.pyc differ
diff --git a/image_processing/__pycache__/utils.cpython-36.pyc b/image_processing/__pycache__/utils.cpython-36.pyc
new file mode 100644
index 0000000..804fc9c
Binary files /dev/null and b/image_processing/__pycache__/utils.cpython-36.pyc differ
diff --git a/image_processing/acquisitionlib.py b/image_processing/acquisitionlib.py
index 7cd6deb..032477a 100755
--- a/image_processing/acquisitionlib.py
+++ b/image_processing/acquisitionlib.py
@@ -68,14 +68,14 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
     nd2_iterator = ND2Reader(nd2file)
     axes = nd2_iterator.axes  # list of char ['t', 'm', 'c', 'z',] 'y' or 'x'
     sizes = nd2_iterator.sizes  # dict with keys in axes list and value is the corresponding dimension
-    print vars(nd2_iterator).keys()
+    print(list(vars(nd2_iterator).keys()))
     bname = os.path.splitext(os.path.basename(nd2_iterator.filename))[0]
-    print "{:<20s}{:<s}".format('bname',bname)
+    print("{:<20s}{:<s}".format('bname',bname))
 
-    print "Axes:", axes
-    print "Sizes:"
+    print("Axes:", axes)
+    print("Sizes:")
     for ax in axes:
-        print "{ax:<5s}{dim:<10d}".format(ax=ax,dim=sizes[ax])
+        print("{ax:<5s}{dim:<10d}".format(ax=ax,dim=sizes[ax]))
 
     # build filtering indexes
     idx = {}
@@ -96,7 +96,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
             tend = nt-1
         tstart = max(tstart,0)
         tend = min(tend,nt-1)
-        print "tstart = {:d}  tend = {:d}".format(tstart,tend)
+        print("tstart = {:d}  tend = {:d}".format(tstart,tend))
         idx['t']=np.arange(tstart,tend+1)
 
     ## fov
@@ -106,9 +106,9 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
         fmtdict['m']=fmt
         if (fovs is None):
             fovs = np.arange(nm)
-        print "FOVs"
+        print("FOVs")
         for i in fovs:
-            print "{:<2s}FOV {:d} selected".format("",i)
+            print("{:<2s}FOV {:d} selected".format("",i))
         idx['m']=fovs
 
     ## colors
@@ -119,7 +119,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
         if (colors is None):
             colors = np.arange(nc)
         for i in colors:
-            print "Channel {:d} selected".format(i)
+            print("Channel {:d} selected".format(i))
         idx['c']=colors
 
     ## cropping
@@ -135,7 +135,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
         ylo=0
         yhi=ny-1
     if not (ylo < yhi):
-        print "Problem with y-cropping: ylo = {:d}    yhi = {:d}".format(ylo,yhi)
+        print("Problem with y-cropping: ylo = {:d}    yhi = {:d}".format(ylo,yhi))
     if np.mod(yhi+1-ylo,2) == 1:
         yhi = yhi - 1
     idx['y']=np.arange(ylo,yhi+1, dtype=np.uint)
@@ -146,12 +146,12 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
         nx = sizes['x']
         xlo,xhi = np.asarray(xcrop,dtype=np.uint)
         xlo = max(xlo,0)
-        xhi = min(xhi,nx-1)
+        xhi = min(xhi,ny-1)
     else:
         xlo=0
         xhi=nx-1
     if not (xlo < xhi):
-        print "Problem with x-cropping: xlo = {:d}    xhi = {:d}".format(xlo,xhi)
+        print("Problem with x-cropping: xlo = {:d}    xhi = {:d}".format(xlo,xhi))
     if np.mod(xhi+1-xlo,2) == 1:
         xhi = xhi - 1
     idx['x']=np.arange(xlo,xhi+1,dtype=np.uint)
@@ -172,14 +172,14 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
 
     if 'm' in axes:
         nd2_iterator.iter_axes='m'
-        print "Starting per-FOV writing"
+        print("Starting per-FOV writing")
         for fov in idx['m']:
-            print "FOV {:d}".format(fov)
+            print("FOV {:d}".format(fov))
             frame = nd2_iterator[fov]
             fov_no = frame.frame_no
             tiff_meta = frame.metadata
             if (fov_no != fov):
-                print "Inconsistency for fov {:d}: fov_no={:d}".format(fov,fov_no)
+                print("Inconsistency for fov {:d}: fov_no={:d}".format(fov,fov_no))
             if 't' in axes:
                 raise ValueError("Time axis handling not implemented yet.")
                 frame.iter_axes='t'
@@ -198,7 +198,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
                 fname = fname.format(fov=fov)
                 fileout = os.path.join(tiffdir,fname+'.tif')
                 ti.imwrite(fileout, img, imagej=True, photometric='minisblack',metadata=tiff_meta)
-                print "{:<20s}{:<s}".format('fileout', fileout)
+                print("{:<20s}{:<s}".format('fileout', fileout))
     else:
         frame = nd2_iterator
         tiff_meta = frame.metadata
@@ -219,7 +219,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
                 # write tiff as a stack
                 fileout = os.path.join(tiffdir,fname.format(t=t) +'.tif')
                 ti.imwrite(fileout, img, imagej=True, photometric='minisblack',metadata=tiff_meta)
-                print "{:<20s}{:<s}".format('fileout', fileout)
+                print("{:<20s}{:<s}".format('fileout', fileout))
             sys.exit()
         else:
             img = np.array(frame[0])
@@ -234,7 +234,7 @@ def process_nd2_tiff(nd2file, tstart=None, tend=None, fovs=None, colors=None, xc
             fname = bname
             fileout = os.path.join(tiffdir,fname+'.tif')
             ti.imwrite(fileout, img, imagej=True, photometric='minisblack',metadata=tiff_meta)
-            print "{:<20s}{:<s}".format('fileout', fileout)
+            print("{:<20s}{:<s}".format('fileout', fileout))
 
     return nd2_iterator.metadata
 
diff --git a/image_processing/collection_cells.py b/image_processing/collection_cells.py
index 41daeb0..9b29d54 100755
--- a/image_processing/collection_cells.py
+++ b/image_processing/collection_cells.py
@@ -9,7 +9,7 @@
 import shutil
 import tifffile as ti
 import cv2
-import cPickle as pkl
+import pickle as pkl
 
 # custom
 from utils import *
@@ -18,18 +18,18 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 #################### function ####################
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['collection']={}
@@ -46,18 +46,18 @@ def check_segmentation(sdir):
     OUTPUT:
       * boolean
     """
-    print "{:<20s}{:<s}".format("segmentation dir:", sdir)
+    print("{:<20s}{:<s}".format("segmentation dir:", sdir))
 
     # check directory existence
     if not os.path.isdir(sdir):
-        print "Segmentation directory doesn\'t exist."
+        print("Segmentation directory doesn\'t exist.")
         return False
 
     # check index
     indexname='index_tiffs.txt'
     pathtoindex=os.path.join(sdir,indexname)
     if not os.path.isfile(pathtoindex):
-        print "Index file doesn\'t exist: {:s}".format(indexname)
+        print("Index file doesn\'t exist: {:s}".format(indexname))
         return False
 
     # check existence of all files in the index
@@ -72,22 +72,22 @@ def check_segmentation(sdir):
 
         # check tiff
         if not check_tiff_file(f):
-            print "Problem with file {:s}".format(f)
+            print("Problem with file {:s}".format(f))
             return False
 
         # check estimator file
         if not check_estimator_file(ef):
-            print "Problem with file {:s}".format(ef)
+            print("Problem with file {:s}".format(ef))
             return False
 
         # check mask file
         if not check_mask_file(mf):
-            print "Problem with file {:s}".format(mf)
+            print("Problem with file {:s}".format(mf))
             return False
 
         # check labels file
         if not check_labels_file(lf):
-            print "Problem with file {:s}".format(lf)
+            print("Problem with file {:s}".format(lf))
             return False
 
     return True
@@ -155,20 +155,20 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
         plt.xticks([]),plt.yticks([])
         fileout = os.path.join(debugdir,bname+".png")
         plt.savefig(fileout, dpi=300, bbox_inches='tight', pad_inches=0)
-        print "{:<20s}{:<s}".format('fileout',fileout)
+        print("{:<20s}{:<s}".format('fileout',fileout))
         plt.close('all')
 
     # write tiff
     fileout = os.path.join(tiff_dir,bname+'.tif')
     ti.imwrite(fileout, subimg, imagej=True, photometric='minisblack')
-    print "{:<20s}{:<s}".format('fileout',fileout)
+    print("{:<20s}{:<s}".format('fileout',fileout))
 
     # write
     fileout = os.path.join(mask_dir,bname+'.tif')
     #fileout = os.path.join(mask_dir,bname+'.txt')
     #np.savetxt(fileout,submask)
     ti.imwrite(fileout, np.array(255*submask,dtype=np.uint8), imagej=True, photometric='minisblack')
-    print "{:<20s}{:<s}".format('fileout',fileout)
+    print("{:<20s}{:<s}".format('fileout',fileout))
 
     return
 
@@ -177,7 +177,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
     parser = argparse.ArgumentParser(prog="Segmentation tool -- cells.")
     parser.add_argument('-f', '--paramfile',  type=file, required=False, help='Yaml file containing parameters.')
     parser.add_argument('-d', '--outputdir',  type=str, required=False, help='Output directory')
-    #parser.add_argument('--lean',  action='store_true', required=False, help='Do not write crops for collection.')
+    parser.add_argument('--lean',  action='store_true', required=False, help='Do not write crops for collection.')
     parser.add_argument('--debug',  action='store_true', required=False, help='Enable debug mode')
 
     # INITIALIZATION
@@ -194,7 +194,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
     outputdir = os.path.join(outputdir,'cells','collection')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
@@ -215,7 +215,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
 
     # debug or not
     if namespace.debug:
-        print "!! Debug mode !!"
+        print("!! Debug mode !!")
 
     # GET METADATA
     pathtometa = os.path.join(rootdir,'metadata.txt')
@@ -228,7 +228,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
     has_segmentation = check_segmentation(seg_dir)
     if not has_segmentation:
         raise ValueError("Problems with segmentation! Run segmentation again.")
-    print "Segmentation detected"
+    print("Segmentation detected")
 
     # parameters
     mpp = params['px2um']
@@ -249,11 +249,11 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
     tiff_dir = os.path.join(outputdir,'tiffs')
     if not os.path.isdir(tiff_dir):
         os.makedirs(tiff_dir)
-    print "{:<20s}{:<s}".format("tiff_dir", tiff_dir)
+    print("{:<20s}{:<s}".format("tiff_dir", tiff_dir))
     mask_dir = os.path.join(outputdir,'masks')
     if not os.path.isdir(mask_dir):
         os.makedirs(mask_dir)
-    print "{:<20s}{:<s}".format("mask_dir", mask_dir)
+    print("{:<20s}{:<s}".format("mask_dir", mask_dir))
     ## tiff list
     indexname='index_tiffs.txt'
     pathtoindex=os.path.join(seg_dir,indexname)
@@ -262,34 +262,34 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
     for n in range(nfiles):
         # load tiff
         f, ef, mf, lf = index[n]
-        print "Processing file {:d} / {:d}".format(n,nfiles)
+        print("Processing file {:d} / {:d}".format(n,nfiles))
         f = os.path.relpath(os.path.join(seg_dir,f))
-        print f
+        print(f)
         with ti.TiffFile(f) as tif:
             img = tif.asarray()
             meta = tif.imagej_metadata
         if mpp is None:
             try:
                 mpp = float(meta['mpp'])
-                print "mpp = {:.6f}".format(mpp)
-            except KeyError, ValueError:
-                print "Unit is lacking: mpp = 1!"
+                print("mpp = {:.6f}".format(mpp))
+            except KeyError as ValueError:
+                print("Unit is lacking: mpp = 1!")
                 mpp = 1.
         fov = meta['m']
         nchannels, height, width = img.shape
-        print img.dtype
-        print img.shape
+        print(img.dtype)
+        print(img.shape)
 
         # make meshgrid
         Y,X = np.mgrid[0:height,0:width]
 
         # load labels
         lf = os.path.relpath(os.path.join(seg_dir,lf))
-        print lf
+        print(lf)
         labels = ssp.load_npz(lf).todense()
         nlabels = np.max(labels)
         if nlabels == 0:
-            print "No labels detected"
+            print("No labels detected")
             continue
 
         # compute background
@@ -301,7 +301,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
         # iterate over segmented objects
         for n in range(1, nlabels):
             cell = {}
-            print "Getting cell {:d} / {:d} for FOV {:d}".format(n,nlabels,fov)
+            print("Getting cell {:d} / {:d} for FOV {:d}".format(n,nlabels,fov))
 
             # make index
             mask = (labels == n)
@@ -326,7 +326,7 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
             cell['bounding_box_rotated']['ycoord_center']=xymid[1]
             cell['bounding_box_rotated']['width']=w
             cell['bounding_box_rotated']['height']=h
-            cell['bounding_box_rotated']['angle']=angle
+            cell['bounding_box_rotated']['angle']=h
 
             # dimensions
             if w > h:
@@ -359,15 +359,15 @@ def write_crop(img, mask, points, bname, tiff_dir='.', mask_dir='.', pad_x=5, pa
             cells.append(cell)
 
             # write tiff
-            if params['write_cropped']:
+            if not namespace.lean:
                 write_crop(img, mask, points, bname=cell_id, tiff_dir=tiff_dir, mask_dir=mask_dir,debug=namespace.debug, **params['crops'])
 
     ncells = len(cells)
-    print "ncells = {:d} collected".format(ncells)
+    print("ncells = {:d} collected".format(ncells))
 
 # write down the cell dictionary
     celldict = {cell['id']: cell for cell in cells}
     celldict = make_dict_serializable(celldict)
     pathtocells = os.path.join(outputdir, 'collection.js')
     write_dict2json(pathtocells,celldict)
-    print "{:<20s}{:<s}".format('fileout', pathtocells)
+    print("{:<20s}{:<s}".format('fileout', pathtocells))
diff --git a/image_processing/preprocess_images.py b/image_processing/preprocess_images.py
index d52cee1..9c6e5de 100755
--- a/image_processing/preprocess_images.py
+++ b/image_processing/preprocess_images.py
@@ -10,7 +10,7 @@
 #from PIL import Image
 import tifffile as ti
 import cv2
-import cPickle as pkl
+import pickle as pkl
 import scipy.ndimage as simg
 
 # custom
@@ -20,18 +20,18 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 #################### function ####################
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['preprocess_images']={}
@@ -65,15 +65,15 @@ def get_background_checkerboard(img, size=65):
     for j in range(nh):
         h0 = hedges[j]
         h1 = hedges[j+1]
-        print "h0 = {:.1e}    h1 = {:.1e}".format(h0,h1)
+        print("h0 = {:.1e}    h1 = {:.1e}".format(h0,h1))
         for i in range(nw):
             w0 = wedges[i]
             w1 = wedges[i+1]
-            print "w0 = {:.1e}    w1 = {:.1e}".format(w0,w1)
+            print("w0 = {:.1e}    w1 = {:.1e}".format(w0,w1))
             rect = img[h0:h1:,w0:w1]
 
             val = np.median(rect)
-            print "median = {:.1e}".format(val)
+            print("median = {:.1e}".format(val))
             bg[h0:h1,w0:w1] = val
             img_coarse[j,i] = val
 
@@ -154,16 +154,6 @@ def get_background_medianblur(img, size=201):
     myimg = (myimg - imin) / (imax - imin)
     myimg = np.array(255*myimg, dtype=np.uint8)
     bg = cv2.medianBlur(myimg, ksize=size)
-#
-#    size = 2* (size//2) + 1 # I like odd numbers
-#    footprint = np.ones((size,size))
-#    nrange = np.arange(size)
-#    NN,MM = np.meshgrid(nrange,nrange)
-#    n0 = size//2
-#    idx = ((NN-n0)**2+(MM-n0)**2 < (0.5*size)**2)
-#    footprint[~idx] = 0
-#    myimg = simg.median_filter(myimg, footprint=footprint, mode='reflect')
-
 
     bg = cv2.blur(bg,ksize=(size,size)) # maybe not essential, to smooth background values
     bg = np.array(bg, dtype=np.float32)/255
@@ -172,108 +162,6 @@ def get_background_medianblur(img, size=201):
 
     return bg
 
-def get_background_logmean_fft(img, size=201, logtransform=False):
-    shape = img.shape
-    if len(shape) != 2:
-        raise ValueError("img must be a 2D array")
-    dtype = img.dtype
-    norm = get_img_norm(dtype)
-
-    # rescale
-    myimg = np.array(img, dtype=np.float32)
-    imin = np.min(myimg)
-    imax = np.max(myimg)
-    scale = (imax-imin)
-
-    offset = 0.
-    if (imin < 0.):
-        offset = 2*np.abs(imin)
-    elif (imin ==  0.):
-        offset = np.sort(np.unique(myimg))[1]
-    myimg = (myimg+offset)/scale
-
-    # do the fft rescaling
-    h,w = myimg.shape
-    l = size//2
-    #print("l = {:d}".format(l))
-    X = np.zeros((h+2*l,w+2*l))
-    hnew,wnew = X.shape
-    F = np.zeros((hnew,wnew))
-    hrange = np.arange(hnew)
-    wrange = np.arange(wnew)
-
-    kfreq_h = np.fft.fftfreq(hnew) # frequencies
-    kfreq_w = np.fft.fftfreq(wnew) # frequencies
-
-    # pad input image by reflection
-    if (l > 0):
-        X[l:-l,l:-l] = np.copy(myimg)
-        Xpad = X[l:2*l,l:-l]
-        X[:l,l:-l] = Xpad[::-1,:]
-        Xpad = X[-2*l:-l,l:-l]
-        X[-l:,l:-l] = Xpad[::-1,:]
-
-        Xpad = X[:,l:2*l]
-        X[:,:l] = Xpad[:,::-1]
-        Xpad = X[:,-2*l:-l]
-        X[:,-l:] = Xpad[:,::-1]
-    else:
-        X = np.copy(myimg)
-
-    # make footprint for average
-    footprint_func = lambda n,m : (n**2 + m**2 <= l**2) | ((hnew-n)**2 + m**2 <= l**2) | ((hnew-n)**2 + (wnew-m)**2 <= l**2) | (n**2 + (wnew-m)**2 <= l**2)
-    hh,ww = np.meshgrid(hrange,wrange,indexing='ij')
-    F = np.float_(footprint_func(hh,ww))
-    fsum = np.sum(F)  # normalization so that it is a mean
-    F /= fsum
-    #print("fsum = {:.6f}".format(fsum))
-    # important to check than without lognorm, with l=0, it reproduces the input after fft and ifft.
-
-    Xcpy = np.copy(X)
-    # log transform
-    if (logtransform):
-        X = np.log(X)
-
-    # fourier transform
-    Xtilde = np.fft.fftn(X)
-    Ftilde = np.fft.fftn(F)
-    X = np.real(np.fft.ifftn(Xtilde*Ftilde))
-    #print(np.max(X-Xcpy), np.median(X-Xcpy))
-
-    # log detransform
-    if (logtransform):
-        X = np.exp(X)
-    bg = X[l:-l,l:-l]
-
-#    # test
-#    import matplotlib.pyplot as plt
-#    xmin = np.min(myimg)
-#    xmax = np.max(myimg)
-#    fig=plt.figure()
-#    ax = fig.add_subplot(1,2,1)
-#    #im = ax.imshow(F)
-#    im = ax.imshow(Xcpy, vmin=xmin,vmax=xmax)
-#    fig.colorbar(im,ax=ax)
-#
-#    ax = fig.add_subplot(1,2,2)
-#    #ax.imshow(np.real(Ftilde))
-#    im = ax.imshow(X, vmin=xmin,vmax=xmax)
-#
-#    fig.colorbar(im,ax=ax)
-#
-#    fig.savefig('footprint.png')
-#    plt.close()
-#    sys.exit()
-#    # test
-
-    # restore scaling
-    bg = cv2.blur(bg,ksize=(size,size)) # maybe not essential, to smooth background values
-    #bg = np.array(bg, dtype=np.float32)/255
-    bg = scale* bg - offset
-    bg = np.array(bg, dtype=img.dtype)
-    return bg
-
-
 def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg_size=200, debug=False):
     """
     INPUT:
@@ -296,9 +184,8 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
     """
 
     # method for background subtraction
-    #get_background = get_background_medianblur
+    get_background = get_background_medianblur
     #get_background = get_background_checkerboard
-    get_background = lambda arr, size: get_background_logmean_fft(arr, size=size, logtransform=True)
 
     # pre-processing
     bname = os.path.splitext(os.path.basename(tiff_file))[0]
@@ -324,11 +211,11 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
     if invert == None:
         invert = []
     elif invert == 'all':
-        invert = range(nchannels)
+        invert = list(range(nchannels))
 
     # adjust background sliding window size
     bg_size = 2*int(bg_size/2) + 1 # make odd
-    print "bg_size = {:d}".format(bg_size)
+    print("bg_size = {:d}".format(bg_size))
 
     img0 = np.copy(img)
     # process channels
@@ -348,7 +235,7 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
         if bg_subtract:
 #            print "Background subtraction..."
             # method for background subtraction using sliding window
-            img_bg[c] = get_background(arr, bg_size)
+            img_bg[c] = get_background(arr, size=bg_size)
             idx = arr > img_bg[c]
 #            print np.unique(arr)
 #            print np.unique(img_bg[c])
@@ -375,15 +262,15 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
     fname = os.path.basename(tiff_file)
     fileout = os.path.join(outputdir,fname)
     ti.imwrite(fileout, img, imagej=True, photometric='minisblack', metadata=meta)
-    print "{:<20s}{:<s}".format('fileout',fileout)
+    print("{:<20s}{:<s}".format('fileout',fileout))
 
     # debug
     if debug:
-        print "Debug for preprocess images"
+        print("Debug for preprocess images")
         img_bg_post = np.zeros(shape, dtype=dtype)
         for c in range(nchannel):
 #            print "c = {:d}".format(c)
-            img_bg_post[c] = get_background(img[c], bg_size)
+            img_bg_post[c] = get_background(img[c], size=bg_size)
 #            print np.unique(img[c])
 #            print np.unique(img_bg_post[c])
         debugdir = os.path.join(outputdir,'debug')
@@ -423,7 +310,7 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
         fileout = os.path.join(debugdir,fname + '.png')
         gs.tight_layout(fig,w_pad=0)
         plt.savefig(fileout,dpi=300)
-        print "{:<20s}{:<s}".format('debug file', fileout)
+        print("{:<20s}{:<s}".format('debug file', fileout))
         plt.close('all')
 
     # exit
@@ -450,7 +337,7 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
         outputdir = os.path.relpath(outputdir, os.getcwd())
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
@@ -482,14 +369,13 @@ def preprocess_image(tiff_file, outputdir='.', invert=None, bg_subtract=True, bg
     # copy metadata
     tiffdir=os.path.dirname(os.path.relpath(tiff_files[0], rootdir))
     metadata=os.path.join(tiffdir,'metadata.txt')
-    if os.path.isfile(metadata):
-        dest = os.path.join(outputdir,os.path.basename(metadata))
-        if (os.path.realpath(metadata) != os.path.realpath(dest)):
-            shutil.copy(metadata,dest)
+    dest = os.path.join(outputdir,os.path.basename(metadata))
+    if (os.path.realpath(metadata) != os.path.realpath(dest)):
+        shutil.copy(metadata,dest)
 
     # debug or not
     if namespace.debug:
-        print "!! Debug mode !!"
+        print("!! Debug mode !!")
 
     params=allparams['preprocess_images']
 
diff --git a/image_processing/process_collection.py b/image_processing/process_collection.py
index 2acd4ac..e56a936 100755
--- a/image_processing/process_collection.py
+++ b/image_processing/process_collection.py
@@ -23,12 +23,12 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 # matplotlib controls
 plt.rcParams['svg.fonttype'] = 'none'  # to embed fonts in output ('path' is to convert as text as paths)
@@ -39,7 +39,7 @@
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['process_collection']={}
@@ -67,7 +67,7 @@ def default_parameters():
 
     cells = load_json2dict(cellfile)
     ncells = len(cells)
-    print "ncells = {:d}".format(ncells)
+    print("ncells = {:d}".format(ncells))
 
     # output directory
     outputdir = namespace.outputdir
@@ -77,7 +77,7 @@ def default_parameters():
         outputdir = os.path.relpath(outputdir, os.getcwd())
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
@@ -98,12 +98,12 @@ def default_parameters():
     params = allparams['process_collection']
     new_attributes = params['new_attributes']
     mpp = params['mpp']
-    keys = cells.keys()
+    keys = list(cells.keys())
     keys_sel = []
     for n in range(ncells):
         key = keys[n]
         cell = cells[key]
-        print "Cell {:d} / {:d} : {:s}".format(n,ncells,cell['id'])
+        print("Cell {:d} / {:d} : {:s}".format(n,ncells,cell['id']))
 
         # filtering operations
 
@@ -119,11 +119,11 @@ def default_parameters():
             if mpp is None:
                 try:
                     mpp = cell['mpp']
-                except KeyError, ValueError:
-                    print "Unit is lacking: mpp = 1! Skipping cell."
+                except KeyError as ValueError:
+                    print("Unit is lacking: mpp = 1! Skipping cell.")
                     continue
             if not 'volume' in cell:
-                print 'Volume must be computed for \'volume_um\''
+                print('Volume must be computed for \'volume_um\'')
                 continue
             cell['volume_um']  = cell['volume']*mpp*mpp*mpp
 
@@ -135,11 +135,11 @@ def default_parameters():
     # WRITE PROCESSED CELL FILE
     cells_new = {key: cells[key] for key in keys_sel}
     ncells_new = len(cells_new)
-    print "ncells_new = {:d}".format(ncells_new)
+    print("ncells_new = {:d}".format(ncells_new))
 
     cells_new = make_dict_serializable(cells_new)
     bname = os.path.splitext(os.path.basename(cellfile))[0]
     fname = "{}_processed".format(bname)
     pathtocells = os.path.join(outputdir, fname + '.js')
     write_dict2json(pathtocells,cells_new)
-    print "{:<20s}{:<s}".format('fileout', pathtocells)
+    print("{:<20s}{:<s}".format('fileout', pathtocells))
diff --git a/image_processing/process_nd2.py b/image_processing/process_nd2.py
index a021945..aba8057 100755
--- a/image_processing/process_nd2.py
+++ b/image_processing/process_nd2.py
@@ -16,18 +16,18 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 #################### function ####################
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['process_nd2']={}
@@ -47,7 +47,7 @@ def default_parameters():
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(prog="Read and handle nd2 files.")
     parser.add_argument('ND2',  type=str, help='ND2 file to open.')
-    parser.add_argument('-f', '--paramfile',  type=file, required=False, help='Yaml file containing parameters.')
+    parser.add_argument('-f', '--paramfile',  type=open, required=False, help='Yaml file containing parameters.')
     parser.add_argument('-d', '--outputdir',  type=str, required=False, help='Output directory')
     parser.add_argument('--debug',  action='store_true', required=False, help='Enable debug mode')
 
@@ -63,7 +63,7 @@ def default_parameters():
     outputdir = os.path.join(outputdir,'TIFFS')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print(("{:<20s}{:<s}".format("outputdir", outputdir)))
 
     # parameter file
     if namespace.paramfile is None:
@@ -91,4 +91,3 @@ def default_parameters():
 #    with open(fileout,'w') as fout:
 #        yaml.safe_dump(metainfo, fout, encoding=('utf-8'), default_flow_style=False, allow_unicode=False)
 #
-
diff --git a/image_processing/segmentation_cells.py b/image_processing/segmentation_cells.py
index c74884b..567194b 100755
--- a/image_processing/segmentation_cells.py
+++ b/image_processing/segmentation_cells.py
@@ -10,7 +10,7 @@
 #from PIL import Image
 import tifffile as ti
 import cv2
-import cPickle as pkl
+import pickle as pkl
 
 # custom
 from utils import *
@@ -19,18 +19,18 @@
 # yaml formats
 npfloat_representer = lambda dumper,value: dumper.represent_float(float(value))
 nparray_representer = lambda dumper,value: dumper.represent_list(value.tolist())
-float_representer = lambda dumper,value: dumper.represent_scalar(u'tag:yaml.org,2002:float', "{:<.8e}".format(value))
+float_representer = lambda dumper,value: dumper.represent_scalar('tag:yaml.org,2002:float', "{:<.8e}".format(value))
 unicode_representer = lambda dumper,value: dumper.represent_unicode(value.encode('utf-8'))
 yaml.add_representer(float,float_representer)
 yaml.add_representer(np.float_,npfloat_representer)
 yaml.add_representer(np.ndarray,nparray_representer)
-yaml.add_representer(unicode,unicode_representer)
+yaml.add_representer(str,unicode_representer)
 
 #################### function ####################
 def default_parameters():
     """Generate a default parameter dictionary."""
 
-    print "Loading default parameters"
+    print("Loading default parameters")
     params={}
     # filtering - select only a subset
     params['segmentation']={}
@@ -157,14 +157,14 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
     # rescale dynamic range linearly (important for OTSU)
     amin = np.min(img)
     amax = np.max(img)
-    print "amin = {:.1g}    amax = {:.1g}".format(amin,amax)
+    print("amin = {:.1g}    amax = {:.1g}".format(amin,amax))
     img = (np.float_(img) - amin)/(amax-amin)
 
     ## thresholding to binary mask
     norm8 = float(2**8-1)
     norm16 = float(2**16-1)
     if threshold is None:
-        print "OTSU thresholding"
+        print("OTSU thresholding")
         # convert to 8-bit image (Open CV requirement for OTSU)
         img = np.array(255*img,np.uint8)
         img8 = np.copy(img)
@@ -179,13 +179,13 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
         thres16 = np.uint16(thres1*norm16)
         ret = max((thres1 - amin)/(amax-amin),0) # value in rescaled DNR
         ret = np.uint8(255*ret) # uint8
-        print "thres1 = {:.1g}    threshold_rescaled_DNR_uint8 = {:d}".format(thres1,ret)
+        print("thres1 = {:.1g}    threshold_rescaled_DNR_uint8 = {:d}".format(thres1,ret))
         # convert to 8-bit image (Open CV requirement for OTSU)
         img = np.array(255*img,np.uint8)
         img8 = np.copy(img)
         # input threshold
         ret1,img = cv2.threshold(img,ret,255,cv2.THRESH_BINARY)
-    print "thres1 = {:.1g}    thres8 = {:d}    thres16 = {:d}".format(thres1, thres8, thres16)
+    print("thres1 = {:.1g}    thres8 = {:d}    thres16 = {:d}".format(thres1, thres8, thres16))
     img_bin = np.copy(img)
 
     ## opening/closing operations
@@ -199,7 +199,7 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
 
     ## find connected components
     ncomp, labels = cv2.connectedComponents(img)
-    print "Found {:d} objects".format(ncomp)
+    print("Found {:d} objects".format(ncomp))
 
     ## compute the bounding box for the identified labels
     #for n in range(ncomp):
@@ -278,7 +278,7 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
         eimg[idx] = scores[n]
     nz = np.sum(eimg > 0.)
     ntot = len(np.ravel(eimg))
-    print "nz = {:d} / {:d}    sparcity index = {:.2e}".format(nz, ntot, float(nz)/float(ntot))
+    print("nz = {:d} / {:d}    sparcity index = {:.2e}".format(nz, ntot, float(nz)/float(ntot)))
     efname = bname
     #efile = os.path.join(outputdir,efname+'.txt')
     #efile = os.path.join(outputdir,efname+'.pkl')
@@ -287,7 +287,7 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
         #np.savetxt(fout, eimg)
         #pkl.dump(eimg,fout)
         ssp.save_npz(efile, ssp.coo_matrix(eimg), compressed=False)
-    print "{:<20s}{:<s}".format('est. file', efile)
+    print("{:<20s}{:<s}".format('est. file', efile))
 
     if debug:
         debugdir = os.path.join(outputdir,'debug')
@@ -359,7 +359,7 @@ def get_estimator_boundingbox(tiff_file, channel=0, outputdir='.', w0=1, w1=100,
         fileout = os.path.join(debugdir,fname + '.png')
         gs.tight_layout(fig,w_pad=0)
         plt.savefig(fileout,dpi=300)
-        print "{:<20s}{:<s}".format('debug file', fileout)
+        print("{:<20s}{:<s}".format('debug file', fileout))
         plt.close('all')
     # from agarpad code: end """
 
@@ -442,7 +442,7 @@ def get_mask(f, ef, threshold=0., outputdir='.', debug=False):
         #np.savetxt(fout, mask)
         #pkl.dump(mask,fout)
         ssp.save_npz(mfile, ssp.coo_matrix(mask), compressed=False)
-    print "{:<20s}{:<s}".format('mask file', mfile)
+    print("{:<20s}{:<s}".format('mask file', mfile))
 
     # debug
     if debug:
@@ -462,7 +462,7 @@ def get_mask(f, ef, threshold=0., outputdir='.', debug=False):
         fileout = os.path.join(debugdir,fname + '.png')
         fig.tight_layout()
         plt.savefig(fileout,dpi=300, bbox_inches='tight', pad_inches=0)
-        print "{:<20s}{:<s}".format('debug file', fileout)
+        print("{:<20s}{:<s}".format('debug file', fileout))
         plt.close('all')
 
     return os.path.realpath(mfile)
@@ -492,7 +492,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         #np.savetxt(fout, labels)
         #pkl.dump(labels,fout)
         ssp.save_npz(lfile, ssp.coo_matrix(labels), compressed=False)
-    print "{:<20s}{:<s}".format('labels file', lfile)
+    print("{:<20s}{:<s}".format('labels file', lfile))
 
     # debug
     if debug:
@@ -514,7 +514,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         fileout = os.path.join(debugdir,fname + '.png')
         fig.tight_layout()
         plt.savefig(fileout,dpi=300, bbox_inches='tight', pad_inches=0)
-        print "{:<20s}{:<s}".format('debug file', fileout)
+        print("{:<20s}{:<s}".format('debug file', fileout))
         plt.close('all')
 
 
@@ -542,7 +542,7 @@ def get_label(f, mf, outputdir='.', debug=False):
     outputdir = os.path.join(outputdir,'cells','segmentation')
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
     # parameter file
     if namespace.paramfile is None:
@@ -571,7 +571,7 @@ def get_label(f, mf, outputdir='.', debug=False):
 
     # debug or not
     if namespace.debug:
-        print "!! Debug mode !!"
+        print("!! Debug mode !!")
 
     params=allparams['segmentation']
     segmentation_method=params['method']
@@ -591,12 +591,12 @@ def get_label(f, mf, outputdir='.', debug=False):
     if ntiffs == 0:
         raise ValueError("Number of tiff files is zero!")
     else:
-        print "{:<20s}{:<d}".format("ntiffs", ntiffs)
+        print("{:<20s}{:<d}".format("ntiffs", ntiffs))
 
     if not os.path.isfile(pathtoindex):
         index = np.asarray([[os.path.relpath(f,outputdir)] for f in tiff_files])
         write_index(pathtoindex,index)
-        print "{:<20s}{:<s}".format("fileout", pathtoindex)
+        print("{:<20s}{:<s}".format("fileout", pathtoindex))
 
     # move metadata if found
     index = load_index(pathtoindex)
@@ -617,7 +617,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         estimator_dir=os.path.join(outputdir,'estimators')
         if not os.path.isdir(estimator_dir):
             os.makedirs(estimator_dir)
-        print "{:<20s}{:<s}".format("est. dir.", estimator_dir)
+        print("{:<20s}{:<s}".format("est. dir.", estimator_dir))
         est_files = []
         for f in tiff_files:
             ef=get_estimator(f, method=segmentation_method, outputdir=estimator_dir, estimator_params=params['estimator_params'][segmentation_method], channel=params['channel'], debug=namespace.debug)
@@ -626,7 +626,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         index = np.concatenate([index,np.transpose([est_files])], axis=1)
 
         write_index(pathtoindex,index)
-        print "{:<20s}{:<s}".format("fileout", pathtoindex)
+        print("{:<20s}{:<s}".format("fileout", pathtoindex))
 
     # BUILD MASKS
     build_mask=False
@@ -638,7 +638,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         mask_dir=os.path.join(outputdir,'masks')
         if not os.path.isdir(mask_dir):
             os.makedirs(mask_dir)
-        print "{:<20s}{:<s}".format("mask dir.", mask_dir)
+        print("{:<20s}{:<s}".format("mask dir.", mask_dir))
         mask_files = []
         nfile = len(index)
         for n in range(nfile):
@@ -652,7 +652,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         index = np.concatenate([index,np.transpose([mask_files])], axis=1)
 
         write_index(pathtoindex,index)
-        print "{:<20s}{:<s}".format("fileout", pathtoindex)
+        print("{:<20s}{:<s}".format("fileout", pathtoindex))
 
     # BUILD LABELS
     build_label=False
@@ -664,7 +664,7 @@ def get_label(f, mf, outputdir='.', debug=False):
         label_dir=os.path.join(outputdir,'labels')
         if not os.path.isdir(label_dir):
             os.makedirs(label_dir)
-        print "{:<20s}{:<s}".format("label dir.", label_dir)
+        print("{:<20s}{:<s}".format("label dir.", label_dir))
         label_files = []
         nfile = len(index)
         for n in range(nfile):
@@ -679,4 +679,4 @@ def get_label(f, mf, outputdir='.', debug=False):
         index = np.concatenate([index,np.transpose([label_files])], axis=1)
 
         write_index(pathtoindex,index)
-        print "{:<20s}{:<s}".format("fileout", pathtoindex)
+        print("{:<20s}{:<s}".format("fileout", pathtoindex))
diff --git a/image_processing/utils.py b/image_processing/utils.py
index 0f2aacd..b2cfedf 100755
--- a/image_processing/utils.py
+++ b/image_processing/utils.py
@@ -105,7 +105,7 @@ def load_json2dict(pathtojson):
     return mydict
 
 def make_dict_serializable(mydict):
-    for k, v in mydict.iteritems():
+    for k, v in mydict.items():
         if isinstance(v, dict):
             make_dict_serializable(v)
         else:
@@ -140,8 +140,8 @@ def get_tiff2ndarray(tiff_file,channel=0,metadata=False, normalize=True):
         with ti.TiffFile(tiff_file) as tif:
             img = tif.asarray()
             meta = tif.imagej_metadata
-    except Exception, e:
-        print e
+    except Exception as e:
+        print(e)
         raise ValueError("Opening tiff with tiff failed.")
 
     naxis = len(img.shape)
@@ -238,7 +238,7 @@ def get_OTSU(tiff_files, outputdir='.', write=True,NBYTESMAX=1000000000, nfiles=
         header="channel, float, uint8, uint16"
         fileout = os.path.join(outputdir, "thresholds_otsu.txt")
         np.savetxt(fileout,thresholds, header=header,fmt='%-8d%-8.4f%-8d%-8d')
-        print "{:<20s}{:<s}".format("fileout",fileout)
+        print("{:<20s}{:<s}".format("fileout",fileout))
 
     # return
     return thresholds
@@ -263,7 +263,7 @@ def get_OTSU(tiff_files, outputdir='.', write=True,NBYTESMAX=1000000000, nfiles=
         outputdir = os.path.relpath(outputdir, os.getcwd())
     if not os.path.isdir(outputdir):
         os.makedirs(outputdir)
-    print "{:<20s}{:<s}".format("outputdir", outputdir)
+    print("{:<20s}{:<s}".format("outputdir", outputdir))
 
 #    # parameter file
 #    if namespace.paramfile is None:
@@ -287,10 +287,10 @@ def get_OTSU(tiff_files, outputdir='.', write=True,NBYTESMAX=1000000000, nfiles=
         if test:
             tiff_files.append(f)
     ntiffs = len(tiff_files)
-    print "ntiffs = {:d}".format(ntiffs)
+    print("ntiffs = {:d}".format(ntiffs))
 
     # OTSU thresholds
     if namespace.otsu and (ntiffs > 0):
-        print "Computing OTSU thresholds..."
+        print("Computing OTSU thresholds...")
         thresholds = get_OTSU(tiff_files, outputdir=outputdir)
 
diff --git a/roles/params.yml b/roles/params.yml
new file mode 100755
index 0000000..92021fd
--- /dev/null
+++ b/roles/params.yml
@@ -0,0 +1,56 @@
+process_nd2:
+  colors:
+  fovs:
+  xcrop: [100,1000]
+  ycrop: [100,1000]
+
+preprocess_images:
+  invert: [2]
+  bg_subtract: True
+  bg_size: 200
+
+# um/px: 0.11
+segmentation:
+  channel: 2
+  method: bounding_box
+  estimator_params:
+    bounding_box:
+      acut: 0.8
+      aratio_min: 2.5
+      aratio_max: 100
+      h0: 20
+      h1: 200
+      w0: 8.5
+      w1: 10.5
+      border_pad: 5
+      threshold:
+  mask_params:
+    threshold: 0.95
+
+collection:
+  # normally the px2um is saved in the metadata.txt file
+  px2um: null
+  crops:
+    pad_x: 5
+    pad_y: 5
+
+dimensions:
+  bins:
+    - 32
+    - 16
+    - 32
+    - 32
+  units_dx:
+    - 2
+    - 0.1
+    - 2
+    - 2
+  mode: um
+
+dimensions_other:
+  bins:
+    - 32
+    - 32
+  units_dx:
+    - 2
+    - 0.1