diff --git a/.gitignore b/.gitignore
index f6b0aa8..3370aba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -89,3 +89,15 @@ ENV/
 .ropeproject
 
 .DS_Store
+
+# imagingdata
+data/
+
+# output
+output/
+
+# jupyter-notebook
+ipynb/
+
+# ~
+*~
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
index 0856c7c..ee67c66 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -13,6 +13,10 @@ RUN conda install matplotlib==1.5.1 \
     pypng==0.0.18 mahotas==1.4.1  opencv-python==3.2.0.7 \
     git+https://github.com/jfrelinger/cython-munkres-wrapper \
     jupyter
+RUN pip install numba notebook==5.4.1
+RUN pip install fast-histogram
+
+
 
 EXPOSE 8888
 WORKDIR /home
diff --git a/celltk/caller.py b/celltk/caller.py
index 2c8b979..b1d4bbf 100644
--- a/celltk/caller.py
+++ b/celltk/caller.py
@@ -1,4 +1,3 @@
-
 import argparse
 from os.path import join, isdir, exists
 from glob import glob
@@ -7,7 +6,6 @@
 import yaml
 import multiprocessing
 from utils.file_io import make_dirs
-print 'test celltk'
 import sys
 
 logger = logging.getLogger(__name__)
@@ -83,7 +81,6 @@ def run_operation(output_dir, operation):
     functions, params, images, labels, output = parse_operation(operation)
     inputs = prepare_path_list(images, output_dir)
     logger.info(inputs)
-    
     inputs_labels = prepare_path_list(labels, output_dir)
     output = join(output_dir, output) if output else output_dir
     caller = _retrieve_caller_based_on_function(functions[0])
diff --git a/celltk/labeledarray/labeledarray/labeledarray.py b/celltk/labeledarray/labeledarray/labeledarray.py
index cb92c95..0ebf573 100644
--- a/celltk/labeledarray/labeledarray/labeledarray.py
+++ b/celltk/labeledarray/labeledarray/labeledarray.py
@@ -83,7 +83,7 @@ def _label2idx(self, item):
         if boolarr.all():
             return (slice(None, None, None), ) + (slice(None, None, None), ) * (self.ndim - 1)
         minidx = min(tidx) if min(tidx) > 0 else None
-        maxidx = max(tidx) if max(tidx) < self.shape[0] - 1 else None
+        maxidx = max(tidx)+1 if max(tidx)+1 < self.shape[0] else None
         if boolarr.sum() > 1:
             return (slice(minidx, maxidx, None), ) + (slice(None, None, None), ) * (self.ndim - 1)
 
diff --git a/celltk/postprocess.py b/celltk/postprocess.py
index b314c71..b4eee73 100644
--- a/celltk/postprocess.py
+++ b/celltk/postprocess.py
@@ -85,7 +85,7 @@ def main():
     parser.add_argument("-l", "--labels", help="labels", nargs="+")
     parser.add_argument("-o", "--output", help="output directory", type=str, default='temp')
     parser.add_argument("-f", "--functions", help="functions", nargs="+")
-    parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+    parser.add_argument('-p', '--param', nargs='+', help='parameters', action='append')
     args = parser.parse_args()
 
     if args.functions is None:
diff --git a/celltk/preprocess.py b/celltk/preprocess.py
index 92efc19..d1a862a 100644
--- a/celltk/preprocess.py
+++ b/celltk/preprocess.py
@@ -37,7 +37,9 @@ def main():
     parser.add_argument("-i", "--input", help="images", nargs="*")
     parser.add_argument("-o", "--output", help="output directory", type=str, default='temp')
     parser.add_argument("-f", "--functions", help="functions", nargs="*")
-    parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+    parser.add_argument('-p', '--param', nargs='+', help='parameters', action='append')
+    # parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+
     args = parser.parse_args()
 
     if args.functions is None:
diff --git a/celltk/preprocess_operation.py b/celltk/preprocess_operation.py
index 0ea4c78..94bd244 100644
--- a/celltk/preprocess_operation.py
+++ b/celltk/preprocess_operation.py
@@ -16,17 +16,17 @@
 from utils.global_holder import holder
 from utils.mi_align import calc_jitters_multiple, calc_crop_coordinates
 from utils.shading_correction import retrieve_ff_ref
-import cv2
-from utils.background_subtractor import subtract_background_rolling_ball
 from scipy.ndimage.filters import gaussian_filter
 
 logger = logging.getLogger(__name__)
 np.random.seed(0)
 
+
 def gaussian_blur(img, SIGMA=3):
     img = gaussian_filter(img, sigma=SIGMA)
     return img
 
+
 def gaussian_laplace(img, SIGMA=2.5, NEG=False):
     if NEG:
         img = -calc_lapgauss(img, SIGMA)
@@ -50,14 +50,15 @@ def background_subtraction_wavelet_hazen(img, THRES=100, ITER=5, WLEVEL=6, OFFSE
     img = img - back
     return convert_positive(img, OFFSET)
 
-def background_subtraction_rolling_ball_hazen(img, RADIUS=100, SIGMA=3, OFFSET=50):
+
+def rolling_ball(img, RADIUS=100, SIGMA=3, OFFSET=50):
     """Rolling ball background subtraction.
     """
     back = rolling_ball_subtraction_hazen(img.astype(np.float), RADIUS)
     img = img - back
-    #return img
     return convert_positive(img, OFFSET)
 
+
 def n4_illum_correction(img, RATIO=1.5, FILTERINGSIZE=50):
     """
     Implementation of the N4 bias field correction algorithm.
@@ -109,33 +110,6 @@ def align(img, CROP=0.05):
         return img[jt[0]:jt[1], jt[2]:jt[3], :]
 
 
-def align2(img, CROP=0.05):
-    """
-    CROP (float): crop images beforehand. When set to 0.05, 5% of each edges are cropped.
-    """
-    if not hasattr(holder, "align"):
-        if isinstance(holder.inputs[0], list) or isinstance(holder.inputs[0], tuple):
-            inputs = [i[0] for i in holder.inputs]
-        else:
-            inputs = holder.inputs
-
-        img0 = imread(inputs[0])
-
-        (ch, cw) = [int(CROP * i) for i in img0.shape]
-        ch = None if ch == 0 else ch
-        cw = None if cw == 0 else cw
-
-        jitters = calc_jitters_multiple(inputs, ch, cw)
-        holder.align = calc_crop_coordinates(jitters, img0.shape)
-        logger.debug('holder.align set to {0}'.format(holder.align))
-    jt = holder.align[holder.frame]
-    logger.debug('Jitter: {0}'.format(jt))
-    if img.ndim == 2:
-        return img[jt[0]:jt[1], jt[2]:jt[3]]
-    if img.ndim == 3:
-        return img[jt[0]:jt[1], jt[2]:jt[3], :]
-
-
 def flatfield_references(img, ff_paths=['Pos0/img00.tif', 'Pos1/img01.tif'], exp_corr=False):
     """
     Use empty images for background subtraction and illumination bias correction.
@@ -227,6 +201,19 @@ def correct_shade(img, ref, darkref, ch):
     img = correct_shade(img, ref, darkref, ch)
     return img
 
+def correct_uneven_illumination(img, bkgimg='Pos1/min_stack.tif'):
+    '''
+    correct_uneven_illumination
+    bkgimg is background image such as blank image or minimum projection result.
+    '''
+    img = img.astype(np.float)
+    bg = gaussian_blur(imread(bkgimg), 3)
+    bg = bg.astype(np.float)
+    d0 = img - bg
+    d0[d0 < 0] = 0
+    m_bg = 1 / (bg / bg.max()) # nega image of background
+    img = d0 * m_bg
+    return img
 
 def background_subtraction_wavelet(img, level=7, OFFSET=10):
     '''
@@ -262,11 +249,28 @@ def stitch_images(img, POINTS=[(0,0),(0,0),(0,0),(0,0)]):
     img = np_arithmetic(img, 'max')
     return img
 
-def rolling_ball(img, RADIUS=30):
-    
-    img = (img/256).astype('uint8')
-    img = subtract_background_rolling_ball(img, RADIUS, light_background=False,\
-                                     use_paraboloid=False, do_presmooth=True, create_background=False)
-    return img
-
-
+def bleedthourh_correction(img, BT=0.0):
+    img0 = img[:, :, 0].astype(np.float)
+    img1 = img[:, :, 1].astype(np.float)
+    img0 = img0 - img1 * BT
+    return img0
+
+def deep_unet(img, weight_path, region=1):
+    """ Generates a probability map of cells using the UNet algorithm. 
+
+    Args:
+        img (numpy.ndarray): image that will be segmented 
+        weight_path (string): path to weights file in .hdf5 format, can either bea local path or url. 
+        region (int): determines if each image will be saved individually (region =1) or 
+            setting region as None will save a stack of float32 images.
+    Returns:
+        pimg (numpy.ndarray): probability map image 
+        
+    """
+    from utils.unet_predict import predict
+    from utils.file_io import LocalPath
+    from utils.global_holder import holder
+    with LocalPath(weight_path) as wpath:
+        pimg = predict(holder.path, wpath)
+    pimg = np.moveaxis(pimg, 0, -1)
+    return pimg[:, :, region]
diff --git a/celltk/segment.py b/celltk/segment.py
index 4583f3b..fa5ae28 100644
--- a/celltk/segment.py
+++ b/celltk/segment.py
@@ -60,11 +60,10 @@ def main():
     parser.add_argument("-o", "--output", help="output directory",
                         type=str, default='temp')
     parser.add_argument("-f", "--functions", help="functions", nargs="*", default=None)
-    parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+    parser.add_argument('-p', '--param', nargs='+', help='parameters', action='append')
     args = parser.parse_args()
 
     params = ParamParser(args.param).run()
-
     if args.functions is None:
         print help(segment_operation)
         return
diff --git a/celltk/segment_operation.py b/celltk/segment_operation.py
index e5b0ea2..7a67fc8 100644
--- a/celltk/segment_operation.py
+++ b/celltk/segment_operation.py
@@ -44,7 +44,7 @@ def adaptive_thres_two(img, FIL1=10, FIL2=100, R1=100, R2=100):
     """
     bw = adaptive_thresh(img, R=R1, FILTERINGSIZE=FIL1)
     foreground = adaptive_thresh(img, R=R2, FILTERINGSIZE=FIL2)
-    bw[-foreground] = 0
+    bw[~foreground] = 0
     return label(bw)
 
 
@@ -54,7 +54,7 @@ def adaptive_thres_otsu(img, FIL1=4, R1=1):
     """
     bw = adaptive_thresh(img, R1, FIL1)
     foreground = global_otsu(img) > 0
-    bw[-foreground] = 0
+    bw[~foreground] = 0
     return label(bw)
 
 
@@ -67,7 +67,7 @@ def watershed_labels(labels, REG=10):
 
 def lap_peak_local(img, separation=10, percentile=64, min_sigma=2, max_sigma=5, num_sigma=10):
     sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
-    gl_images = [-gaussian_laplace(img, s) * s ** 2 for s in sigma_list]
+    gl_images = [~gaussian_laplace(img, s) * s ** 2 for s in sigma_list]
     image_cube = np.dstack(gl_images)
     max_image = np.max(image_cube, axis=2)
     coords = grey_dilation(max_image, separation=separation, percentile=percentile)
@@ -81,20 +81,3 @@ def mark_pos(im, coords):
     return label(binary_dilation(bw, np.ones((3, 3))))
 
 
-def deepcell(img, model_path, weight_path, padding=30, rad=[10, 30]):
-    """
-    model_path and weight_path can be either local path or url.
-    """
-    from utils.tf_deepcell.predict import predict
-    from segment import clean_labels
-    from subdetect_operation import propagate_multisnakes
-    from utils.file_io import LocalPath
-    from utils.global_holder import holder
-    with LocalPath(model_path) as mpath, LocalPath(weight_path) as wpath:
-        pimg = predict(holder.path, mpath, wpath)
-    cell = pimg[1] > pimg[0] * 100
-    cell[pimg[1] < pimg[2] * 100] = False
-    cell = np.pad(cell, padding, 'constant')
-    print cell.shape, img.shape
-    cimg = propagate_multisnakes(label(cell), img, NITER=2, lambda2=30)
-    return clean_labels(cimg, rad=rad)
diff --git a/celltk/subdetect.py b/celltk/subdetect.py
index ccfc529..3577cd4 100644
--- a/celltk/subdetect.py
+++ b/celltk/subdetect.py
@@ -43,7 +43,7 @@ def main():
     parser.add_argument("-o", "--output", help="output directory",
                         type=str, default='temp')
     parser.add_argument("-f", "--functions", help="functions", nargs="+")
-    parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+    parser.add_argument('-p', '--param', nargs='+', help='parameters', action='append')
     args = parser.parse_args()
 
     if args.functions is None:
diff --git a/celltk/subdetect_operation.py b/celltk/subdetect_operation.py
index f6f80bb..7d078df 100644
--- a/celltk/subdetect_operation.py
+++ b/celltk/subdetect_operation.py
@@ -8,6 +8,8 @@
 import numpy as np
 from scipy.ndimage import morphology
 from skimage.morphology import remove_small_objects
+from utils.labels_handling import convert_labels
+from utils.subdetect_utils import label_high_pass, label_nearest
 
 np.random.seed(0)
 
@@ -130,7 +132,7 @@ def morphological(labels, func='grey_opening', size=3, iterations=1):
 def watershed_divide(labels, regmax=10, min_size=100):
     """
     divide objects in labels with watershed segmentation.
-        regmax: 
+        regmax:
         min_size: objects smaller than this size will not be divided.
     """
     from utils.subdetect_utils import watershed_labels
@@ -141,3 +143,68 @@ def watershed_divide(labels, regmax=10, min_size=100):
     ws_large += labels.max()
     ws_large[ws_large == labels.max()] = 0
     return labels + ws_large
+
+
+def cytoplasm_levelset(labels, img, niter=20, dt=-0.5, thres=0.5):
+    """ Segment cytoplasm from supplied nuclear labels and probability map 
+        Expand using level sets method from nuclei to membrane.
+        Uses an implementation of Level set method. See: 
+        https://wiseodd.github.io/techblog/2016/11/05/levelset-method/
+    
+    Args:
+        labels (numpy.ndarray): nuclear mask labels
+        img (numpy.ndarray): probability map 
+        niter (int): step size to expand mask, number of iterations to run the levelset algorithm
+        dt (float): negative values for porgation, positive values for shrinking 
+        thres (float): threshold of probability value to extend the nuclear mask to
+    
+    Returns:
+        cytolabels (numpy.ndarray): cytoplasm mask labels  
+
+    """
+    from skimage.morphology import closing, disk, remove_small_holes
+    from utils.dlevel_set import dlevel_set
+    phi = labels.copy()
+    phi[labels == 0] = 1
+    phi[labels > 0] = -1
+
+    outlines = img.copy()
+    outlines = -outlines
+    outlines = outlines - outlines.min()
+    outlines = outlines/outlines.max()
+
+    mask = outlines < thres
+    phi = dlevel_set(phi, outlines, niter=niter, dt=dt, mask=mask)
+
+    cytolabels = label(remove_small_holes(label(phi < 0)))
+    cytolabels = closing(cytolabels, disk(3))
+
+    temp = cytolabels.copy()
+    temp[labels == 0] = 0
+    cytolabels = convert_labels(temp, labels, cytolabels)
+    return cytolabels
+
+def segment_bacteria(nuc, img, slen=3, SIGMA=0.5,THRES=100, CLOSE=20, THRESCHANGE=1000, MINAREA=5):
+   """ Segment bacteria and assign to closest nucleus
+
+    Args:
+        nuc (numpy.ndarray): nuclear mask labels
+        img (numpy.ndarray): image in bacterial channel
+        slen (int): Size of Gaussian kernel
+        SIGMA (float): Standard deviation for Gaussian kernel
+        THRES (int): Threshold pixel intensity fo real signal 
+        CLOSE (int): Radius for disk used to return morphological closing of the image (dilation followed by erosion to remove dark spots and connect bright cracks)
+        THRESCHANGE (int): argument unnecessary? 
+        MINAREA (int): minimum area in pixels for a bacterium 
+
+    Returns:
+        labels (numpy.ndarray[np.uint16]): bacterial mask labels  
+
+    """
+   labels = label_high_pass(img, slen=slen, SIGMA=SIGMA, THRES=THRES, CLOSE=3)
+   if labels.any():
+           labels, comb, nuc_prop, nuc_loc = label_nearest(img, labels, nuc)
+   from skimage.morphology import remove_small_objects
+   labels = remove_small_objects(labels, MINAREA)
+   return labels.astype(np.uint16)
+
diff --git a/celltk/track.py b/celltk/track.py
index 2a6d4e5..ed25d07 100644
--- a/celltk/track.py
+++ b/celltk/track.py
@@ -57,7 +57,7 @@ def main():
     parser.add_argument("-o", "--output", help="output directory",
                         type=str, default='temp')
     parser.add_argument("-f", "--functions", help="functions", nargs="+")
-    parser.add_argument("-p", "--param", nargs="*", help="parameters", default=[])
+    parser.add_argument('-p', '--param', nargs='+', help='parameters', action='append')
     args = parser.parse_args()
 
     if args.functions is None:
diff --git a/celltk/track_operation.py b/celltk/track_operation.py
index c2fd43a..6487251 100644
--- a/celltk/track_operation.py
+++ b/celltk/track_operation.py
@@ -120,7 +120,7 @@ def run_lap(img0, img1, labels0, labels1, DISPLACEMENT=30, MASSTHRES=0.2):
 
 
 def track_neck_cut(img0, img1, labels0, labels1, DISPLACEMENT=10, MASSTHRES=0.2,
-                   EDGELEN=5, THRES_ANGLE=180, WSLIMIT=False, SMALL_RAD=3, CANDS_LIMIT=300):
+                   EDGELEN=5, THRES_ANGLE=180, WSLIMIT=False, SMALL_RAD=None, CANDS_LIMIT=300):
     """
     Adaptive segmentation by using tracking informaiton.
     Separate two objects by making a cut at the deflection. For each points on the outline,
diff --git a/celltk/utils/_model_builder.py b/celltk/utils/_model_builder.py
new file mode 100644
index 0000000..bc04a9d
--- /dev/null
+++ b/celltk/utils/_model_builder.py
@@ -0,0 +1,99 @@
+import keras.layers
+import keras.models
+import tensorflow as tf
+
+CONST_DO_RATE = 0.5
+
+option_dict_conv = {"activation": "relu", "border_mode": "same"}
+option_dict_bn = {"mode": 0, "momentum" : 0.9}
+
+
+# returns a core model from input to 64 channels of the same size
+def get_core(dim1, dim2, dim3):
+
+    # assume dim1 x dim2 image with dim3 color channels
+
+    x = keras.layers.Input(shape=(dim1, dim2, dim3))
+
+    a = keras.layers.Convolution2D(64, 3, 3, **option_dict_conv)(x)
+    a = keras.layers.BatchNormalization(**option_dict_bn)(a)
+
+    a = keras.layers.Convolution2D(64, 3, 3, **option_dict_conv)(a)
+    a = keras.layers.BatchNormalization(**option_dict_bn)(a)
+
+
+    y = keras.layers.MaxPooling2D()(a)
+
+    b = keras.layers.Convolution2D(128, 3, 3, **option_dict_conv)(y)
+    b = keras.layers.BatchNormalization(**option_dict_bn)(b)
+
+    b = keras.layers.Convolution2D(128, 3, 3, **option_dict_conv)(b)
+    b = keras.layers.BatchNormalization(**option_dict_bn)(b)
+
+
+    y = keras.layers.MaxPooling2D()(b)
+
+    c = keras.layers.Convolution2D(256, 3, 3, **option_dict_conv)(y)
+    c = keras.layers.BatchNormalization(**option_dict_bn)(c)
+
+    c = keras.layers.Convolution2D(256, 3, 3, **option_dict_conv)(c)
+    c = keras.layers.BatchNormalization(**option_dict_bn)(c)
+
+
+    y = keras.layers.MaxPooling2D()(c)
+
+    d = keras.layers.Convolution2D(512, 3, 3, **option_dict_conv)(y)
+    d = keras.layers.BatchNormalization(**option_dict_bn)(d)
+
+    d = keras.layers.Convolution2D(512, 3, 3, **option_dict_conv)(d)
+    d = keras.layers.BatchNormalization(**option_dict_bn)(d)
+
+
+    # UP
+
+    d = keras.layers.UpSampling2D()(d)
+
+    y = keras.layers.merge([d, c], concat_axis=3, mode="concat")
+
+    e = keras.layers.Convolution2D(256, 3, 3, **option_dict_conv)(y)
+    e = keras.layers.BatchNormalization(**option_dict_bn)(e)
+
+    e = keras.layers.Convolution2D(256, 3, 3, **option_dict_conv)(e)
+    e = keras.layers.BatchNormalization(**option_dict_bn)(e)
+
+    e = keras.layers.UpSampling2D()(e)
+
+
+    y = keras.layers.merge([e, b], concat_axis=3, mode="concat")
+
+    f = keras.layers.Convolution2D(128, 3, 3, **option_dict_conv)(y)
+    f = keras.layers.BatchNormalization(**option_dict_bn)(f)
+
+    f = keras.layers.Convolution2D(128, 3, 3, **option_dict_conv)(f)
+    f = keras.layers.BatchNormalization(**option_dict_bn)(f)
+
+    f = keras.layers.UpSampling2D()(f)
+
+
+    y = keras.layers.merge([f, a], concat_axis=3, mode="concat")
+
+    y = keras.layers.Convolution2D(64, 3, 3, **option_dict_conv)(y)
+    y = keras.layers.BatchNormalization(**option_dict_bn)(y)
+
+    y = keras.layers.Convolution2D(64, 3, 3, **option_dict_conv)(y)
+    y = keras.layers.BatchNormalization(**option_dict_bn)(y)
+
+    return [x, y]
+
+def get_model(dim1, dim2, dim3, activation="softmax"):
+
+    [x, y] = get_core(dim1, dim2, dim3)
+
+    y = keras.layers.Convolution2D(3, 1, 1, **option_dict_conv)(y)
+
+    if activation is not None:
+        y = keras.layers.Activation(activation)(y)
+
+    model = keras.models.Model(x, y)
+
+    return model
\ No newline at end of file
diff --git a/celltk/utils/_mutinfo.py b/celltk/utils/_mutinfo.py
new file mode 100644
index 0000000..afd6952
--- /dev/null
+++ b/celltk/utils/_mutinfo.py
@@ -0,0 +1,121 @@
+"""
+Replace np histogram to fast-histogram.
+"""
+
+# Copyright (C) 2013 Oskar Maier
+# 
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+# 
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# author Oskar Maier
+# version r0.1.0
+# since 2013-07-09
+# status Release
+
+# build-in modules
+
+# third-party modules
+import numpy
+from fast_histogram import histogram1d, histogram2d
+# code
+def mutual_information(i1, i2, bins=256):
+    r"""
+    Computes the mutual information (MI) (a measure of entropy) between two images.
+
+    MI is not real metric, but a symmetric and nonnegative similarity measures that
+    takes high values for similar images. Negative values are also possible.
+    
+    Intuitively, mutual information measures the information that ``i1`` and ``i2`` share: it
+    measures how much knowing one of these variables reduces uncertainty about the other.
+    
+    The Entropy is defined as:
+    
+    .. math::
+    
+        H(X) = - \sum_i p(g_i) * ln(p(g_i)
+
+    with :math:`p(g_i)` being the intensity probability of the images grey value :math:`g_i`.
+    
+    Assuming two images :math:`R` and :math:`T`, the mutual information is then computed by comparing the
+    images entropy values (i.e. a measure how well-structured the common histogram is).
+    The distance metric is then calculated as follows:
+    
+    .. math::
+    
+        MI(R,T) = H(R) + H(T) - H(R,T) = H(R) - H(R|T) = H(T) - H(T|R)
+    
+    A maximization of the mutual information is equal to a minimization of the joint
+    entropy.
+    
+    Parameters
+    ----------
+    i1 : array_like
+        The first image.
+    i2 : array_like
+        The second image.
+    bins : integer
+        The number of histogram bins (squared for the joined histogram).
+    
+    Returns
+    -------
+    mutual_information : float
+        The mutual information distance value between the supplied images.
+    
+    Raises
+    ------
+    ArgumentError
+        If the supplied arrays are of different shape.
+    """
+    # pre-process function arguments
+    i1 = numpy.asarray(i1)
+    i2 = numpy.asarray(i2)
+    
+    # # validate function arguments
+    # if not i1.shape == i2.shape:
+    #     raise ArgumentError('the two supplied array-like sequences i1 and i2 must be of the same shape')
+    
+    # compute i1 and i2 histogram range
+    i1_range = __range(i1, bins)
+    i2_range = __range(i2, bins)
+    
+    # compute joined and separated normed histograms
+    i1i2_hist = histogram2d(i1.flatten(), i2.flatten(), bins=bins, range=[i1_range, i2_range]) # Note: histogram2d does not flatten array on its own
+    i1_hist = histogram1d(i1, bins=bins, range=i1_range)
+    i2_hist = histogram1d(i2, bins=bins, range=i2_range)
+    
+    # compute joined and separated entropy
+    i1i2_entropy = __entropy(i1i2_hist)
+    i1_entropy = __entropy(i1_hist)
+    i2_entropy = __entropy(i2_hist)
+    
+    # compute and return the mutual information distance
+    return i1_entropy + i2_entropy - i1i2_entropy
+
+def __range(a, bins):
+    '''Compute the histogram range of the values in the array a according to
+    scipy.stats.histogram.'''
+    a = numpy.asarray(a)
+    a_max = a.max()
+    a_min = a.min()
+    s = 0.5 * (a_max - a_min) / float(bins - 1)
+    return (a_min - s, a_max + s)
+ 
+def __entropy(data):
+    '''Compute entropy of the flattened data set (e.g. a density distribution).'''
+    # normalize and convert to float
+    data = data/float(numpy.sum(data))
+    # for each grey-value g with a probability p(g) = 0, the entropy is defined as 0, therefore we remove these values and also flatten the histogram
+    data = data[numpy.nonzero(data)]
+    # compute entropy
+    return -1. * numpy.sum(data * numpy.log2(data))
+    
\ No newline at end of file
diff --git a/celltk/utils/background_subtractor.py b/celltk/utils/background_subtractor.py
deleted file mode 100644
index 9b6a591..0000000
--- a/celltk/utils/background_subtractor.py
+++ /dev/null
@@ -1,474 +0,0 @@
-import cv2
-import numpy as np
-
-"""
-Fully Ported to Python from ImageJ's Background Subtractor.
-Only works for 8-bit greyscale images currently.
-Based on the concept of the rolling ball algorithm described
-in Stanley Sternberg's article,
-"Biomedical Image Processing", IEEE Computer, January 1983.
-Imagine that the 2D grayscale image has a third (height) dimension by the image
-value at every point in the image, creating a surface. A ball of given radius
-is rolled over the bottom side of this surface; the hull of the volume
-reachable by the ball is the background.
-http://rsbweb.nih.gov/ij/developer/source/ij/plugin/filter/BackgroundSubtracter.java.html
-"""
-
-
-def subtract_background_rolling_ball(img, radius, light_background=True,
-                                     use_paraboloid=False, do_presmooth=True,
-                                     create_background=False):
-    """
-    Calculates and subtracts or creates background from image.
-    Arguments:
-    :param img - uint8 np array representing image
-    :param radius - Radius of the rolling ball creating the background (actually a
-                      paraboloid of rotation with the same curvature)
-    :param light_background - Whether the image has a light background.
-    :param do_presmooth - Whether the image should be smoothened (3x3 mean) before creating
-                      the background. With smoothing, the background will not necessarily
-                      be below the image data.
-    :param create_background - Whether to create a background, not to subtract it.
-    :param use_paraboloid - Whether to use the "sliding paraboloid" algorithm.
-    :return img - uint8 np array representing background subtracted image
-    """
-    bs = BackgroundSubtract()
-    return bs.rolling_ball_background(img, radius, light_background, use_paraboloid, do_presmooth, create_background)
-
-
-class BackgroundSubtract:
-    X_DIRECTION = 0
-    Y_DIRECTION = 1
-    DIAGONAL_1A = 2
-    DIAGONAL_1B = 3
-    DIAGONAL_2A = 4
-    DIAGONAL_2B = 5
-
-    def __init__(self):
-        self.width = 0
-        self.height = 0
-
-        self.s_width = 0
-        self.s_height = 0
-
-    def rolling_ball_background(self, img, radius, light_background=True,
-                                use_paraboloid=False, do_presmooth=True,
-                                create_background=False):
-        """
-        Calculates and subtracts background from array.
-        Arguments:
-        :param img - uint8 np array representing image
-        :param radius - Radius of the rolling ball creating the background (actually a
-                          paraboloid of rotation with the same curvature)
-        :param light_background - Whether the image has a light background.
-        :param do_presmooth - Whether the image should be smoothened (3x3 mean) before creating
-                          the background. With smoothing, the background will not necessarily
-                          be below the image data.
-        :param create_background - Whether to create a background, not to subtract it.
-        :param use_paraboloid - Whether to use the "sliding paraboloid" algorithm.
-        :return img - uint8 np array representing background subtracted image
-        """
-        self.height, self.width = img.shape
-        self.s_height, self.s_width = img.shape
-
-        _img = img.copy()
-        if do_presmooth:
-            _img = self._smooth(_img)
-
-        _img = _img.reshape(self.height * self.width)
-
-        invert = False
-        if light_background:
-            invert = True
-
-        ball = None
-        if not use_paraboloid:
-            ball = RollingBall(radius)
-
-        float_img = _img.astype('float64')
-        if use_paraboloid:
-            float_img = self._sliding_paraboloid_float_background(float_img, radius, invert)
-        else:
-            float_img = self._rolling_ball_float_background(float_img, invert, ball)
-
-        if create_background:
-            return float_img.astype('uint8').reshape((self.height, self.width))
-
-        offset = 255.5 if invert else 0.5
-        for p in range(0, self.width*self.height):
-            value = (_img[p]&0xff) - float_img[p] + offset
-            value = max((value, 0))
-            value = min((value, 255))
-            img[int(p / self.width), int(p % self.width)] = value
-        return img
-
-    def _smooth(self, img, window=3):
-        """
-        Applies a 3x3 mean filter to specified array.
-        """
-        kernel = np.ones((window, window), np.float64) / (window*window)
-        img = cv2.filter2D(img, -1, kernel)
-        return img
-
-    def _rolling_ball_float_background(self, float_img, invert, ball):
-        shrink = ball.shrink_factor > 1
-        if invert:
-            float_img = 255 - float_img
-
-        small_img = self._shrink_image(float_img, ball.shrink_factor) if shrink else float_img
-        self._roll_ball(ball, small_img)
-
-        if shrink:
-            float_img = self._enlarge_image(small_img, float_img, ball.shrink_factor)
-
-        if invert:
-            float_img = 255 - float_img
-        return float_img
-
-    def _roll_ball(self, ball, float_img):
-        height, width = self.s_height, self.s_width
-        z_ball = ball.data
-        ball_width = ball.width
-        radius = int(ball_width / 2)
-        cache = [0] * (width * ball_width)
-
-        for y in range(-radius, height + radius):
-            next_line_to_write = (y + radius) % ball_width
-            next_line_to_read = y + radius
-            if next_line_to_read < height:
-                src = next_line_to_read * width
-                dest = next_line_to_write * width
-                cache[dest:dest+width] = float_img[src:src+width]
-                float_img[src:src+width] = float("-inf")
-
-            y0 = max((0, y - radius))
-            y_ball0 = y0 - y + radius
-            y_end = y + radius
-            if y_end >= height:
-                y_end = height - 1
-            for x in range(-radius, width+radius):
-                z = float("inf")
-                x0 = max((0, x - radius))
-                x_ball0 = x0 - x + radius
-                x_end = x + radius
-                if x_end >= width:
-                    x_end = width - 1
-
-                y_ball = y_ball0
-                for yp in range(y0, y_end + 1):
-                    cache_pointer = (yp % ball_width) * width + x0
-                    bp = x_ball0 + y_ball * ball_width
-                    for xp in range(x0, x_end + 1):
-                        z_reduced = cache[cache_pointer] - z_ball[bp]
-                        if z > z_reduced:
-                            z = z_reduced
-                        cache_pointer += 1
-                        bp += 1
-                    y_ball += 1
-
-                y_ball = y_ball0
-                for yp in range(y0, y_end + 1):
-                    p = x0 + yp * width
-                    bp = x_ball0 + y_ball * ball_width
-                    for xp in range(x0, x_end + 1):
-                        z_min = z + z_ball[bp]
-                        if float_img[p] < z_min:
-                            float_img[p] = z_min
-                        p += 1
-                        bp += 1
-                    y_ball += 1
-
-    def _shrink_image(self, img, shrink_factor):
-        height, width = self.height, self.width
-
-        self.s_height, self.s_width = int(height / shrink_factor), int(width / shrink_factor)
-
-        img_copy = img.reshape((height, width)).copy()
-        small_img = np.ones((self.s_height, self.s_width), np.float64)
-
-        for y in range(0, self.s_height):
-            for x in range(0, self.s_width):
-                x_mask_min = shrink_factor * x
-                y_mask_min = shrink_factor * y
-                min_value = img_copy[y_mask_min:y_mask_min + shrink_factor,
-                            x_mask_min:x_mask_min + shrink_factor].min()
-                small_img[y, x] = min_value
-        return small_img.reshape(self.s_height * self.s_width)
-
-    def _enlarge_image(self, small_img, float_img, shrink_factor):
-        height, width = self.height, self.width
-        s_height, s_width = self.s_height, self.s_width
-
-        x_s_indices, x_weigths = self._make_interpolation_arrays(width, s_width, shrink_factor)
-        y_s_indices, y_weights = self._make_interpolation_arrays(height, s_height, shrink_factor)
-        line0 = [0.0] * width
-        line1 = [0.0] * width
-        for x in range(0, width):
-            line1[x] = small_img[x_s_indices[x]] * x_weigths[x] + \
-                       small_img[x_s_indices[x] + 1] * (1.0 - x_weigths[x])
-        y_s_line0 = -1
-        for y in range(0, height):
-            if y_s_line0 < y_s_indices[y]:
-                line0, line1 = line1, line0
-                y_s_line0 += 1
-                s_y_ptr = int((y_s_indices[y] + 1) * s_width)
-                for x in range(0, width):
-                    line1[x] = small_img[s_y_ptr + x_s_indices[x]] * x_weigths[x] + \
-                               small_img[s_y_ptr + x_s_indices[x] + 1] * (1.0 - x_weigths[x])
-            weight = y_weights[y]
-            p = y * width
-            for x in range(0, width):
-                float_img[p] = line0[x] * weight + line1[x] * (1.0 - weight)
-
-                p += 1
-        return float_img
-
-    def _make_interpolation_arrays(self, length, s_length, shrink_factor):
-        s_indices = [0] * length
-        weights = [0.0] * length
-        for i in range(0, length):
-            s_idx = int((i - shrink_factor / 2) / shrink_factor)
-            if s_idx >= s_length - 1:
-                s_idx = s_length - 2
-            s_indices[i] = s_idx
-            distance = (i + 0.5) / shrink_factor - (s_idx + 0.5)
-            weights[i] = 1.0 - distance
-        return s_indices, weights
-
-    def _sliding_paraboloid_float_background(self, float_img, radius, invert):
-        height, width = self.height, self.width
-        cache = [0.0] * max((height, width))
-        next_point = [0] * max((height, width))
-        coeff2 = np.float64(0.5) / radius
-        coeff2_diag = np.float64(1.0) / radius
-
-        if invert:
-            float_img = 255 - float_img
-
-        self._correct_corners(float_img, coeff2, cache, next_point)
-        self._filter1d(float_img, self.X_DIRECTION, coeff2, cache, next_point)
-        self._filter1d(float_img, self.Y_DIRECTION, coeff2, cache, next_point)
-        self._filter1d(float_img, self.X_DIRECTION, coeff2, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_1A, coeff2_diag, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_1B, coeff2_diag, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_2A, coeff2_diag, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_2B, coeff2_diag, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_1A, coeff2_diag, cache, next_point)
-        self._filter1d(float_img, self.DIAGONAL_1B, coeff2_diag, cache, next_point)
-
-        if invert:
-            float_img = 255 - float_img
-
-        return float_img
-
-    def _correct_corners(self, float_img, coeff2, cache, next_point):
-        height, width = self.height, self.width
-        corners = [0] * 4
-        corrected_edges = [0, 0]
-        corrected_edges = self._line_slide_parabola(float_img, 0, 1, width, coeff2, cache, next_point, corrected_edges)
-        corners[0] = corrected_edges[0]
-        corners[1] = corrected_edges[1]
-        corrected_edges = self._line_slide_parabola(float_img, (height - 1) * width, 1, width, coeff2, cache, next_point, corrected_edges)
-        corners[2] = corrected_edges[0]
-        corners[3] = corrected_edges[1]
-        corrected_edges = self._line_slide_parabola(float_img, 0, width, height, coeff2, cache, next_point, corrected_edges)
-        corners[0] += corrected_edges[0]
-        corners[2] += corrected_edges[1]
-        corrected_edges = self._line_slide_parabola(float_img, width - 1, width, height, coeff2, cache, next_point, corrected_edges)
-        corners[1] += corrected_edges[0]
-        corners[3] += corrected_edges[1]
-        diag_length = min((width, height))
-        coeff2_diag = 2 * coeff2
-        corrected_edges = self._line_slide_parabola(float_img, 0, 1 + width, diag_length, coeff2_diag, cache, next_point, corrected_edges)
-        corners[0] += corrected_edges[0]
-        corrected_edges = self._line_slide_parabola(float_img, width - 1, -1 + width, diag_length, coeff2_diag, cache, next_point, corrected_edges)
-        corners[1] += corrected_edges[0]
-        corrected_edges = self._line_slide_parabola(float_img, (height - 1) * width, 1 - width, diag_length, coeff2_diag, cache, next_point, corrected_edges)
-        corners[2] += corrected_edges[0]
-        corrected_edges = self._line_slide_parabola(float_img, width * height - 1, -1 - width, diag_length, coeff2_diag, cache, next_point, corrected_edges)
-        corners[3] += corrected_edges[0]
-
-        float_img[0] = min((float_img[0], corners[0] / 3))
-        float_img[width-1] = min((float_img[width-1], corners[1] / 3))
-        float_img[(height-1)*width] = min((float_img[(height-1)*width], corners[2] / 3))
-        float_img[width*height-1] = min((float_img[width*height-1], corners[3] / 3))
-
-    def _line_slide_parabola(self, float_img, start, inc, length, coeff2, cache, next_point, corrected_edges):
-        min_value = float("inf")
-        last_point = 0
-        first_corner, last_corner = length - 1, 0
-        v_prev1, v_prev2 = 0., 0.
-        curvature_test = 1.999 * coeff2
-
-        p = start
-        for i in range(length):
-            v = float_img[p]
-            cache[i] = v
-            min_value = min((min_value, v))
-            if i >= 2 and v_prev1 + v_prev1- v_prev2 - v < curvature_test:
-                next_point[last_point] = i - 1
-                last_point = i - 1
-            v_prev2 = v_prev1
-            v_prev1 = v
-
-            p += inc
-
-        next_point[last_point] = length - 1
-        next_point[length - 1] = float("inf")
-
-        i1 = 0
-        while i1 < length - 1:
-            v1 = cache[i1]
-            min_slope = float("inf")
-            i2 = 0
-            search_to = length
-            recalculate_limit_now = 0
-
-            j = next_point[i1]
-            while j < search_to:
-                v2 = cache[j]
-                slope = (v2 - v1) / (j - i1) + coeff2 * (j - i1)
-                if slope < min_slope:
-                    min_slope = slope
-                    i2 = j
-                    recalculate_limit_now = -3
-                if recalculate_limit_now == 0:
-                    b = 0.5 * min_slope / coeff2
-                    max_search = i1 + int(b + np.sqrt(b*b + (v1 - min_value) / coeff2) + 1)
-                    if 0 < max_search < search_to:
-                        search_to = max_search
-
-                j = next_point[j]
-                recalculate_limit_now += 1
-
-            if i1 == 0:
-                first_corner = i2
-            if i2 == length - 1:
-                last_corner = i1
-            p = start + (i1 + 1) * inc
-            for j in range(i1 + 1, i2):
-                float_img[p] = v1 + (j - i1) * (min_slope - (j - i1) * coeff2)
-
-                p += inc
-            i1 = i2
-        if corrected_edges is not None:
-            if 4 * first_corner >= length:
-                first_corner = 0
-            if 4 * (length - 1 - last_corner) >= length:
-                last_corner = length - 1
-            v1 = cache[first_corner]
-            v2 = cache[last_corner]
-            slope = (v2 - v1) / (last_corner - first_corner)
-            value0 = v1 - slope * first_corner
-            coeff6 = 0
-            mid = 0.5 * (last_corner + first_corner)
-            for i in range(int((length + 2) / 3), int(2 * length / 3) + 1):
-                dx = (i - mid) * 2 / (last_corner - first_corner)
-                poly6 = dx*dx*dx*dx*dx*dx - 1
-                if cache[i] < value0 + slope*i + coeff6*poly6:
-                    coeff6 = -(value0 + slope*i - cache[i]) / poly6
-            dx = (first_corner - mid) * 2.0 / (last_corner - first_corner)
-            corrected_edges[0] = value0 + coeff6*(dx*dx*dx*dx*dx*dx - 1.0) + coeff2*first_corner*first_corner
-            dx = (last_corner-mid)*2.0/(last_corner-first_corner)
-            corrected_edges[1] = value0 + (length-1)*slope + coeff6*(dx*dx*dx*dx*dx*dx - 1.0) + \
-                                 coeff2*(length-1-last_corner)*(length-1-last_corner)
-        return corrected_edges
-
-    def _filter1d(self, float_img, direction, coeff2, cache, next_point):
-        height, width = self.height, self.width
-        start_line = 0
-        n_lines = 0
-        line_inc = 0
-        point_inc = 0
-        length = 0
-
-        if direction == self.X_DIRECTION:
-            n_lines = height
-            line_inc = width
-            point_inc = 1
-            length = width
-        elif direction == self.Y_DIRECTION:
-            n_lines = width
-            line_inc = 1
-            point_inc = width
-            length = height
-        elif direction == self.DIAGONAL_1A:
-            n_lines = width - 2
-            line_inc = 1
-            point_inc = width + 1
-        elif direction == self.DIAGONAL_1B:
-            start_line = 1
-            n_lines = height - 2
-            line_inc = width
-            point_inc = width + 1
-        elif direction == self.DIAGONAL_2A:
-            start_line = 2
-            n_lines = width
-            line_inc = 1
-            point_inc = width - 1
-        elif direction == self.DIAGONAL_2B:
-            start_line = 0
-            n_lines = height - 2
-            line_inc = width
-            point_inc = width - 1
-        for i in range(start_line, n_lines):
-            start_pixel = i * line_inc
-            if direction == self.DIAGONAL_2B:
-                start_pixel += width - 1
-            if direction == self.DIAGONAL_1A:
-                length = min((height, width-i))
-            elif direction == self.DIAGONAL_1B:
-                length = min((width, height-i))
-            elif direction == self.DIAGONAL_2A:
-                length = min((height, i+1))
-            elif direction == self.DIAGONAL_2B:
-                length = min((width, height-i))
-            self._line_slide_parabola(float_img, start_pixel, point_inc, length, coeff2, cache, next_point, None)
-
-
-class RollingBall:
-    """
-        A rolling ball (or actually a square part thereof)
-        Here it is also determined whether to shrink the image
-    """
-    def __init__(self, radius):
-
-        self.data = []
-        self.width = 0
-
-        if radius <= 10:
-            self.shrink_factor = 1
-            arc_trim_per = 24
-        elif radius <= 30:
-            self.shrink_factor = 2
-            arc_trim_per = 24
-        elif radius <= 100:
-            self.shrink_factor = 4
-            arc_trim_per = 32
-        else:
-            self.shrink_factor = 8
-            arc_trim_per = 40
-        self.build(radius, arc_trim_per)
-
-    def build(self, ball_radius, arc_trim_per):
-        small_ball_radius = ball_radius / self.shrink_factor
-
-        if small_ball_radius < 1:
-            small_ball_radius = 1
-
-        r_square = small_ball_radius * small_ball_radius
-        x_trim = int(arc_trim_per * small_ball_radius / 100)
-        half_width = round(small_ball_radius - x_trim)
-        self.width = int(2 * half_width + 1)
-        self.data = [0] * self.width * self.width
-
-        p = 0
-        for y in range(self.width):
-            for x in range(self.width):
-                x_val = x - half_width
-                y_val = y - half_width
-
-                temp = r_square - x_val * x_val - y_val * y_val
-                self.data[p] = np.sqrt(temp) if temp > 0 else 0
-
-                p += 1
diff --git a/celltk/utils/dlevel_set.py b/celltk/utils/dlevel_set.py
new file mode 100644
index 0000000..b942d05
--- /dev/null
+++ b/celltk/utils/dlevel_set.py
@@ -0,0 +1,49 @@
+import numpy as np
+from subdetect_utils import calc_mask_exclude_overlap
+from filters import label
+
+
+def grad(x):
+    return np.array(np.gradient(x))
+
+
+def norm(x, axis=0):
+    if x.ndim == 3:
+        return np.sqrt(np.sum(np.square(x), axis=axis))
+    else:
+        return np.abs(x)
+
+
+def div(fx, fy):
+    fyy, fyx = grad(fy)
+    fxy, fxx = grad(fx)
+    return fxx + fyy
+
+
+def dot(x, y, axis=0):
+    return np.sum(x * y, axis=axis)
+
+
+def dlevel_set(phi, F,  niter=20, dt=-0.5, mask=None):
+    """
+    Implementation of Level set method. https://wiseodd.github.io/techblog/2016/11/05/levelset-method/
+    It has an extra repulsion and masking processes.
+        niter (int): iteration
+        dt (float): negative values for propagation, positive values for shrinking
+        phi (ndarray): labels, set inner objects as -1 and outside as 1.
+        F (ndarray): img, set boundaries close to 0 and elsewhere close to 1.
+        mask (ndarray[np.bool]): contours do not cross this line.
+    """
+    if mask is None:
+        mask = np.zeros(phi.shape, np.bool)
+    for i in range(niter):
+        dphi = grad(phi)
+        dphi_norm = norm(dphi)
+
+        region = calc_mask_exclude_overlap(label(phi < 0), 2)  # added for absolute repulsion
+        dphi_norm[region] = 0  # added for absolute repulsion
+        dphi_norm[mask] = 0  # added constraints based on the image
+
+        dphi_t = F * dphi_norm
+        phi = phi + dt * dphi_t
+    return phi
\ No newline at end of file
diff --git a/celltk/utils/labels_handling.py b/celltk/utils/labels_handling.py
index 0f80485..38b704b 100644
--- a/celltk/utils/labels_handling.py
+++ b/celltk/utils/labels_handling.py
@@ -17,21 +17,22 @@ def labels_map(lb0, lb1):
 
 
 def convert_labels(lb_ref_to, lb_ref_from, lb_convert):
-    """
-    lb_ref_to: 
-    lb_ref_from: 
-    lb_convert: a labeled image to be converted.
+    """ Maps a reference set of labels onto another set of labels for the same image 
+         Example: Separate segmentation of two objects in an image (nucleus and cytoplasm) leads
+         to different labels and can reconcile the same cell to have the same label for both 
+         objects. 
+    Args:
+        lb_ref_to (np.ndarray): image with object labels to convert to 
+        lb_ref_from (np.ndarray): image with object labels to be converted 
+        lb_convert (np.ndarray): the labeled image to be converted
+
+    Returns:
+        arr (numpy.ndarray): converted labels for objects with reference set of labels 
+
     """
     lbmap_to, lbmap_from = zip(*labels_map(lb_ref_to, lb_ref_from))
     arr = np.zeros(lb_convert.shape, dtype=np.uint16)
     lb = lb_convert.copy()
-
     for n0, n1 in zip(lbmap_to, lbmap_from):
-        arr[lb == n1] = n0
-        lb[lb == n1] = 0
-
-    remained = label(lb)
-    remained = remained + arr.max()
-    remained[remained == arr.max()] = 0
-    arr = arr + remained
-    return arr
\ No newline at end of file
+        arr[lb == n0] = n1
+    return arr
diff --git a/celltk/utils/miopt_align.py b/celltk/utils/miopt_align.py
index 853e566..0b7e01d 100644
--- a/celltk/utils/miopt_align.py
+++ b/celltk/utils/miopt_align.py
@@ -1,5 +1,8 @@
-from medpy.metric.image import mutual_information
+from __future__ import print_function
 import SimpleITK as sitk
+from _mutinfo import mutual_information
+import numpy as np
+
 
 def offset_slice(pixels1, pixels2, i, j):
     '''Return two sliced arrays where the first slice is offset by i,j
@@ -69,12 +72,25 @@ def calc_crop_coordinates(store, shapes):
 
 
 def sitk_translation(img0, img1, off0=0, off1=0):
+    """
+
+
+
+                + off1
+                ^
+                |
+    off0 + <--- + ---> - off0
+                |
+                v
+                - off1
+
+    """
     s0, s1 = sitk.GetImageFromArray(img0), sitk.GetImageFromArray(img1)
     s0, s1 = sitk.Cast(s0, sitk.sitkFloat32), sitk.Cast(s1, sitk.sitkFloat32)
 
     R = sitk.ImageRegistrationMethod()
     R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=250)
-    R.SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200 )
+    R.SetOptimizerAsRegularStepGradientDescent(4.0, .01, 200)
 
     cc = sitk.TranslationTransform(s0.GetDimension())
     cc.SetOffset([-off1, -off0])
@@ -96,4 +112,38 @@ def register_multiseeds(img0, img1, bins=250, initial=(-30, 0, 30)):
             p2, p1 = offset_slice(img1, img0, s0, s1)
             store.append(((s0, s1), mutual_information(p1, p2, bins)))
     store.sort(key=lambda x: x[1])
-    return store[-1]
\ No newline at end of file
+    return store[-1]
+
+
+def _crop_sample_two(num, x, y, patch_h=200, patch_w=200):
+    """
+    """
+    patch_h = patch_h if patch_h % 2 else patch_h + 1
+    patch_w = patch_w if patch_w % 2 else patch_w + 1
+    coords = ([np.random.randint(patch_h, x.shape[0] - patch_h) for i in range(num)], 
+              [np.random.randint(patch_w, x.shape[1] - patch_w) for i in range(num)])
+    h, w = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
+    xstack = np.zeros((num, patch_h, patch_w), np.float32)
+    ystack = xstack.copy()
+    for n, (ch, cw) in enumerate(zip(*coords)):
+        xstack[n, :, :] = x[ch-h:ch+h+1, cw-w:cw+w+1]
+        ystack[n, :, :] = y[ch-h:ch+h+1, cw-w:cw+w+1]
+    return xstack, ystack
+
+
+def register_multiseeds_crop(img0, img1, bins=250, initial=(-30, 0, 30), num=5, patch_size=200):
+    """
+    Align two images and return jitters and the highest mutual information it finds.
+    It will crop an image into a (patch_size x patch_size) randomly and search for jitters
+    until it finds the same values twice or it reaches num attempts.
+    patch_size at least 150 is recommended.
+    """
+    img0s, img1s = _crop_sample_two(num, img0, img1, patch_size, patch_size)
+    st = []
+    for img0, img1 in zip(img0s, img1s):
+        res = register_multiseeds(img0, img1, bins=bins, initial=initial)
+        if [i for i in st if res[0] == i[0]]:
+            return res
+        st.append(res)
+    print('No overlap. It may not be the optimal alignment')
+    return sorted(st, key=lambda x:x[1])[-1]
diff --git a/celltk/utils/parser.py b/celltk/utils/parser.py
index 11115dd..b603dfb 100644
--- a/celltk/utils/parser.py
+++ b/celltk/utils/parser.py
@@ -1,12 +1,60 @@
 import re
 import ast
 import __builtin__
+import ast
+import argparse
+import distutils.util
+
+def split_params(inputs):
+    if "/" not in inputs:
+        return [inputs]
+    store = inputs.split('/')
+    return store
 
 
 class ParamParser(object):
     def __init__(self, param_args):
         self.param_args = param_args
 
+    def run(self):
+        if self.param_args is None:
+            return [{}]
+        parameters = split_params(self.param_args[0])
+        return [self.convert2dict(parameters) for p in parameters]
+
+    def convert2dict(self, param):
+        dictargs={}
+        for p in param:
+            param_kv = p.split('=')
+            dictargs[param_kv[0]] = param_kv[1]
+        for key, value in dictargs.iteritems():
+            value = str(value)
+            if value[0].isdigit():
+                dictargs[key] = ast.literal_eval(value)
+            elif value[0]=='[':
+                try:
+                    dictargs[key] = ast.literal_eval(value)
+                except:
+                    temp = []
+                    for i in value[1:-1].split(','):
+                        if i.isdigit():
+                            temp.append(ast.literal_eval(i))
+                        else:
+                            temp.append(i)
+                    dictargs[key] = temp
+
+            if value == 'True':
+                dictargs[key] = bool(1)
+            elif value == 'False':
+                dictargs[key] = bool(0)
+
+        return dictargs
+        
+
+class ParamParser1(object):
+    def __init__(self, param_args):
+        self.param_args = param_args
+
     def run(self):
         params = self.split_params(self.param_args[:])
         params = self.iter_combine_list(params)
@@ -91,4 +139,4 @@ def parse_image_files(inputs):
         else:
             li.append(element)
     store.append(li)
-    return zip(*store)
\ No newline at end of file
+    return zip(*store)
diff --git a/celltk/utils/preprocess_utils.py b/celltk/utils/preprocess_utils.py
index 691ef2b..d5aadde 100644
--- a/celltk/utils/preprocess_utils.py
+++ b/celltk/utils/preprocess_utils.py
@@ -118,6 +118,7 @@ def resize_img(himg, origshape):
     resized += minh
     return resized
 
+
 def rolling_ball_subtraction_hazen(img, RADIUS=10, SIGMA=3):
     rb = PyRollingBall(ball_radius=RADIUS, smoothing_sigma=SIGMA)
     f = partial(rb.estimateBG)
diff --git a/celltk/utils/shading_correction.py b/celltk/utils/shading_correction.py
index e648df3..3431f5a 100644
--- a/celltk/utils/shading_correction.py
+++ b/celltk/utils/shading_correction.py
@@ -46,7 +46,9 @@ def shading_correction_folder(inputfolder, outputfolder, binning=3, magnificatio
     parentfolder = inputfolder
     for dirname, subdirlist, filelist in os.walk(parentfolder):
         if 'metadata.txt' in filelist:
-            outputdir = join(outputfolder, dirname.split(parentfolder)[-1])
+            sfol_name = dirname.split(parentfolder)[-1]
+            sfol_name = sfol_name if not sfol_name.startswith('/') else sfol_name[1:]
+            outputdir = join(outputfolder, sfol_name)
             if not os.path.exists(outputdir):
                 os.makedirs(outputdir)
             with open(join(dirname, 'metadata.txt')) as mfile:
@@ -55,9 +57,9 @@ def shading_correction_folder(inputfolder, outputfolder, binning=3, magnificatio
             for chnum, ch in enumerate(channels):
                 pathlist = glob(join(dirname, '*channel{0:03d}*'.format(chnum)))
                 for path in pathlist:
-                    if ch == 'PHASE':
-                        img = imread(path)
-                        tiff.imsave(join(outputdir, os.path.basename(path)), img.astype(np.float32))
-                    else:
+                    try:
                         img = correct_shade(imread(path), ref, darkref, ch)
-                        tiff.imsave(join(outputdir, os.path.basename(path)), img.astype(np.float32))
+                    except:
+                        print "ch might not exist as a reference."
+                        img = imread(path)
+                    tiff.imsave(join(outputdir, os.path.basename(path)), img.astype(np.float32))
diff --git a/celltk/utils/subdetect_utils.py b/celltk/utils/subdetect_utils.py
index f4d7ee5..33ac5e2 100644
--- a/celltk/utils/subdetect_utils.py
+++ b/celltk/utils/subdetect_utils.py
@@ -4,6 +4,12 @@
 from skimage.morphology import watershed as skiwatershed
 from skimage.measure import label
 from skimage.feature import peak_local_max
+from skimage.morphology import disk
+from skimage.morphology import closing
+from skimage.measure import label as skim_label
+from scipy.ndimage.filters import gaussian_filter
+from scipy.signal import convolve2d
+from skimage.measure import regionprops
 
 
 def dilate_sitk(labels, RAD):
@@ -134,3 +140,58 @@ def watershed_labels(labels, regmax):
     wlabels[wlabels == labels.max()] = 0
     all_labels = label(labels + wlabels)
     return all_labels
+
+def pairwise_distance(loc1, loc2):
+    xprev = [i[0] for i in loc1]
+    yprev = [i[1] for i in loc1]
+    xcurr = [i[0] for i in loc2]
+    ycurr = [i[1] for i in loc2]
+    xprevTile = np.tile(xprev, (len(xcurr), 1))
+    yprevTile = np.tile(yprev, (len(ycurr), 1))
+    return abs(xprevTile.T - xcurr) + abs(yprevTile.T - ycurr)
+
+def skilabel(bw, conn=2):
+    '''original label might label any objects at top left as 1. To get around this pad it first.'''
+    bw = np.pad(bw, pad_width=1, mode='constant', constant_values=False)
+    label = skim_label(bw, connectivity=conn)
+    label = label[1:-1, 1:-1]
+    return label
+
+def calc_high_pass_kernel(slen, SIGMA):
+    """For Salmonella"""
+    temp = np.zeros((slen, slen))
+    temp[int(slen/2), int(slen/2)] = 1
+    gf = gaussian_filter(temp, SIGMA)
+    norm = np.ones((slen, slen))/(slen**2)
+    return gf - norm
+
+def calc_high_pass(img, slen, SIGMA):
+    """For Salmonella"""
+    kernel = calc_high_pass_kernel(slen, SIGMA)
+    cc = convolve2d(img, kernel, mode='same')
+    return cc
+
+def label_high_pass(img, slen=3, SIGMA=0.5, THRES=50, CLOSE=3):
+    """For Salmonella"""
+    cc = calc_high_pass(img, slen, SIGMA)
+    cc[cc < 0] = 0
+    la = skilabel(cc > THRES, conn=1)
+    la = closing(la, disk(CLOSE))
+    return la
+
+
+def label_nearest(img, label, nuc, DISTHRES=25):
+    """Label objects to the nearest nuc.
+    """
+    nuc_prop = regionprops(nuc, img, cache=False)
+    sal_prop = regionprops(label, img, cache=False)
+    nuc_loc = [i.centroid for i in regionprops(nuc, img, cache=False)]
+    sal_loc = [i.centroid for i in regionprops(label, img, cache=False)]
+    dist = pairwise_distance(nuc_loc, sal_loc)
+    min_dist_arg = np.argmin(dist, axis=0)
+    template = np.zeros(img.shape, np.uint16)
+    for num, (idx, sal) in enumerate(zip(min_dist_arg, sal_prop)):
+        if dist[idx, num] < DISTHRES:
+            template[sal.coords[:, 0], sal.coords[:, 1]] = nuc_prop[idx].label
+    comb = np.max(np.dstack((template, nuc)), axis=2).astype(np.uint16)
+    return template, comb, nuc_prop, nuc_loc
diff --git a/celltk/utils/tf_deepcell/README.md b/celltk/utils/tf_deepcell/README.md
deleted file mode 100644
index a082078..0000000
--- a/celltk/utils/tf_deepcell/README.md
+++ /dev/null
@@ -1,30 +0,0 @@
-# tf_deepcell
-
-#### A simple training and prediction
-```
-python train.py -i data/nuc0.png -l data/labels0.tif -m data/tests_model.py -o output -n 3000 -e 2 -p 61
-python predict.py -i data/nuc1.png -w output/cnn_model_weights.hdf5 -m data/tests_model.py -o output
-```
-
-#### Resume training
-```
-python train.py -i data/nuc0.png -l data/labels0.tif -m data/weights.tests.hdf5
-```
-
-#### Using multiple channels
-```
-python train.py -i data/FIXME -l data/FIXME -m data/tests_model.py
-```
-
-#### Training using multiple images
-```
-python train.py -i data/nuc0.png / data/nuc0.png -l data/labels0.tif / data/labels1.tif -m data/tests_model.py
-```
-
-#### Prediction example
-The following hdf5 was trained with ```-n 1500000 -e 10 -p 61 -b 256``` with GPU.
-```
-python predict.py -i data/nuc1.png -w data/tests_pretrained.hdf5 -m data/tests_model.py -o output
-```
-
-Use tensorflow (1.3.0) and Cuda 8.0
diff --git a/celltk/utils/tf_deepcell/__init__.py b/celltk/utils/tf_deepcell/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/celltk/utils/tf_deepcell/_dilated_pool.py b/celltk/utils/tf_deepcell/_dilated_pool.py
deleted file mode 100644
index 48b8868..0000000
--- a/celltk/utils/tf_deepcell/_dilated_pool.py
+++ /dev/null
@@ -1,73 +0,0 @@
-from tensorflow.python.ops import nn
-from tensorflow import nn
-
-from tensorflow.python.layers.pooling import _Pooling2D
-from tensorflow.python.layers import utils
-try:
-    from tensorflow.python.keras.layers import Layer
-except:
-    from tensorflow.contrib.keras.python.keras.layers import Layer
-
-# from tensorflow.python.keras._impl.keras.utils import conv_utils
-
-
-# class DilatedMaxPool2D(_Pooling2D, Layer):
-#     def __init__(self, pool_size=(2, 2), strides=None, padding='valid',
-#                  data_format='channels_last', name=None,
-#                  dilation_rate=2, **kwargs):
-#         self.dilation_rate = dilation_rate
-#         if strides is None or dilation_rate > 1:
-#             strides = (1, 1)
-#         super(DilatedMaxPool2D, self).__init__(
-#               nn.pool, pool_size=pool_size, strides=strides,
-#               padding=padding, data_format=data_format, name=name, **kwargs)
-
-#     def call(self, inputs):
-#         outputs = self.pool_function(
-#             inputs, window_shape=self.pool_size, pooling_type="MAX",
-#             strides=self.strides, padding=self.padding.upper(),
-#             dilation_rate=(self.dilation_rate, self.dilation_rate),
-#             data_format=utils.convert_data_format(self.data_format, 4))
-#         return outputs
-
-#     def get_config(self):
-#         config = {
-#             'pool_size': self.pool_size,
-#             'padding': self.padding,
-#             'strides': self.strides,
-#             'data_format': self.data_format,
-#             'dilation_rate': self.dilation_rate,
-#         }
-#         base_config = super(DilatedMaxPool2D, self).get_config()
-#         return dict(list(base_config.items()) + list(config.items()))
-
-
-class DilatedMaxPool2D(_Pooling2D, Layer):
-    def __init__(self, pool_size=(2, 2), strides=None, padding='valid',
-                 data_format='channels_last', name=None,
-                 dilation_rate=2, **kwargs):
-        self.dilation_rate = dilation_rate
-        if strides is None or dilation_rate > 1:
-            strides = (1, 1)
-        super(DilatedMaxPool2D, self).__init__(
-              nn.pool, pool_size=pool_size, strides=strides,
-              padding=padding, data_format=data_format, name=name, **kwargs)
-
-    def call(self, inputs):
-        outputs = self.pool_function(
-            inputs, window_shape=self.pool_size, pooling_type="MAX",
-            strides=self.strides, padding=self.padding.upper(),
-            dilation_rate=(self.dilation_rate, self.dilation_rate),
-            data_format=utils.convert_data_format(self.data_format, 4))
-        return outputs
-
-    def get_config(self):
-        config = {
-            'pool_size': self.pool_size,
-            'padding': self.padding,
-            'strides': self.strides,
-            'data_format': self.data_format,
-            'dilation_rate': self.dilation_rate,
-        }
-        base_config = super(DilatedMaxPool2D, self).get_config()
-        return dict(list(base_config.items()) + list(config.items()))
diff --git a/celltk/utils/tf_deepcell/patches.py b/celltk/utils/tf_deepcell/patches.py
deleted file mode 100644
index b17b361..0000000
--- a/celltk/utils/tf_deepcell/patches.py
+++ /dev/null
@@ -1,194 +0,0 @@
-from __future__ import division
-import numpy as np
-import os
-try:
-    from tensorflow.python.keras import backend
-    from tensorflow.python.keras.preprocessing.image import ImageDataGenerator, Iterator, array_to_img
-except:
-    from tensorflow.contrib.keras.python.keras.preprocessing.image import ImageDataGenerator, Iterator, array_to_img
-    from tensorflow.contrib.keras.python.keras import backend as K
-
-
-def _sample_coords_weighted(num, shape, weights):
-    flat_idx = np.arange(shape[0] * shape[1])
-    chosen_flat = np.random.choice(flat_idx, num, p=weights/weights.sum())
-    return np.unravel_index(chosen_flat, shape)
-
-
-def _calc_equal_weights(features):
-    ap = []
-    for i in np.unique(features):
-        ap.append((features == i).sum())
-    frac = np.array([(np.sum(ap) - i)/np.sum(ap) for i in ap])
-    prob_2d = np.zeros(features.shape)
-    for i in np.unique(features):
-        prob_2d[features == i] = frac[i]
-    return prob_2d
-
-
-def pick_coords(num, features, patch_h, patch_w):
-    """
-    features: img with labels
-    """
-    prob_2d = _calc_equal_weights(features.astype(np.uint8))
-    _ph, _pw = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
-    perim = np.zeros(prob_2d.shape, dtype=np.bool)
-    perim[_ph:-_ph, _pw:-_pw] = True
-    prob_2d[~perim] = 0
-    return _sample_coords_weighted(num, features.shape, prob_2d.flatten())
-
-
-def pick_coords_list(num, li_features, patch_h, patch_w):
-    """
-    features: img with labels
-    """
-    from random import shuffle
-    num = int(num/len(li_features))
-    li_coords = []
-    for en, features in enumerate(li_features):
-        coords = pick_coords(num, features, patch_h, patch_w)
-        coords = zip(*(np.ones(num, np.uint8) * en, coords[0], coords[1]))
-        li_coords.extend(coords)
-    shuffle(li_coords)
-    return li_coords
-
-
-
-# def _pick_coords(num, features, patch_h, patch_w):
-#     """
-#     features: img with labels
-#     """
-#     prob_2d = _calc_equal_weights(features.astype(np.uint8))
-#     _ph, _pw = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
-#     perim = np.zeros(prob_2d.shape, dtype=np.bool)
-#     perim[_ph:-_ph, _pw:-_pw] = True
-#     prob_2d[~perim] = 0
-#     return _sample_coords_weighted(num, features.shape, prob_2d.flatten())
-
-
-def extract_patches(num, x, y, patch_h, patch_w):
-    """
-    x: input images
-    y: feature image with labels
-
-    Sample many windows from a large image. It will correct for labels unbalance.
-    (If there are less labels for cell boundaries, it increases the sampling probability)
-    """
-    coords = pick_coords(num, y, patch_h, patch_w)
-    h, w = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
-    xstack = np.zeros((num, patch_h, patch_w, x.shape[-1]), np.float32)
-    ystack = np.zeros(num)
-    for n, (ch, cw) in enumerate(zip(*coords)):
-        xstack[n, :, :] = x[ch-h:ch+h+1, cw-w:cw+w+1]
-        ystack[n] = y[ch, cw]
-    return xstack, ystack
-
-
-def _extract_patches(x, y, coords, patch_h, patch_w):
-    """
-    x: input images x.shape = (hpix, wpix, ch)
-    y: feature image with labels
-
-    Sample many windows from a large image. It will correct for labels unbalance.
-    (If there are less labels for cell boundaries, it increases the sampling probability)
-    """
-    h, w = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
-    xstack = np.zeros((len(coords), patch_h, patch_w, x.shape[-1]), np.float32)
-    ystack = np.zeros(len(coords))
-    for n, (ch, cw) in enumerate(coords):
-        xstack[n, :, :, :] = x[ch-h:ch+h+1, cw-w:cw+w+1]
-        ystack[n] = y[ch, cw]
-    return xstack, ystack
-
-
-def extract_patch_list(lix, liy, ecoords, patch_h, patch_w):
-    h, w = int(np.floor(patch_h/2)), int(np.floor(patch_w/2))
-    xstack = np.zeros((len(ecoords), patch_h, patch_w, lix[0].shape[-1]), np.float32)
-    ystack = np.zeros(len(ecoords))
-    for n, (cn, ch, cw) in enumerate(ecoords):
-        xstack[n, :, :, :] = lix[cn][ch-h:ch+h+1, cw-w:cw+w+1]
-        ystack[n] = liy[cn][ch, cw]
-    return xstack, ystack
-
-
-class PatchDataGenerator(ImageDataGenerator):
-    def flow(self, x, y, coords, patch_h, patch_w, batch_size=32, shuffle=True, seed=None,
-             save_to_dir=None, save_prefix='', save_format='png'):
-        return CropIterator(
-            x, y, self, coords=coords, patch_h=patch_h, patch_w=patch_w,
-            batch_size=batch_size, shuffle=shuffle, seed=seed,
-            data_format=self.data_format,
-            save_to_dir=save_to_dir, save_prefix=save_prefix, save_format=save_format)
-
-
-class CropIterator(Iterator):
-    def __init__(self, x, y, image_data_generator, coords, patch_h, patch_w,
-                 batch_size=32, shuffle=False, seed=None, func_patch=_extract_patches,
-                 data_format=None, save_to_dir=None, save_prefix='', save_format='png'):
-        self.coords = coords
-        self.patch_h = patch_h
-        self.patch_w = patch_w
-        if isinstance(x, list):
-            self._x, self._y = x[:], y[:]
-            chnum = self._x[0].shape[-1]
-        else:
-            self._x, self._y = x.copy(), y.copy()
-            self._x = self._x[0, :, :, :]
-            chnum = self._x.shape[-1]
-        self.x = np.asarray(np.zeros((1, patch_h, patch_w, chnum)), dtype=K.floatx())
-        if y is not None:
-            self.y = np.zeros(1)
-        else:
-            self.y = None
-        self.n = len(self.coords)
-        self.func_patch = func_patch
-
-        self.image_data_generator = image_data_generator
-        self.data_format = data_format
-        self.save_to_dir = save_to_dir
-        self.save_prefix = save_prefix
-        self.save_format = save_format
-        super(CropIterator, self).__init__(len(self.coords), batch_size, shuffle, seed)
-
-    def _get_batches_of_transformed_samples(self, index_array):
-        batch_x = np.zeros(tuple([len(index_array)] + list(self.x.shape)[1:]), dtype=K.floatx())
-        batch_coords = [self.coords[i] for i in index_array]
-        x, y = self.func_patch(self._x, self._y, batch_coords, self.patch_h, self.patch_w)
-        self.x = x
-        self.y = y
-        index_array = np.arange(len(self.y))
-
-        for i, j in enumerate(index_array):
-            x = self.x[j]
-            x = self.image_data_generator.random_transform(x.astype(K.floatx()))
-            x = self.image_data_generator.standardize(x)
-            batch_x[i] = x
-        if self.save_to_dir:
-            for i, j in enumerate(index_array):
-                img = array_to_img(batch_x[i], self.data_format, scale=True)
-                fname = '{prefix}_{index}_{hash}.{format}'.format(prefix=self.save_prefix, index=j,
-                hash=np.random.randint(1e4), format=self.save_format)
-            img.save(os.path.join(self.save_to_dir, fname))
-        if self.y is None:
-            return batch_x
-        batch_y = self.y[index_array]
-        return batch_x, batch_y
-
-    def next(self):
-        with self.lock:
-            index_array = next(self.index_generator)
-        if isinstance(index_array, tuple):
-            index_array = index_array[0]
-        return self._get_batches_of_transformed_samples(index_array)
-
-
-class PatchDataGeneratorList(ImageDataGenerator):
-    def flow(self, x, y, coords, patch_h, patch_w, batch_size=32, shuffle=True, seed=None,
-             save_to_dir=None, save_prefix='', save_format='png'):
-        x = [i[0] for i in x]
-        return CropIterator(
-            x, y, self, coords=coords, patch_h=patch_h, patch_w=patch_w,
-            batch_size=batch_size, shuffle=shuffle, seed=seed,
-            data_format=self.data_format, func_patch=extract_patch_list,
-            save_to_dir=save_to_dir, save_prefix=save_prefix, save_format=save_format)
-
diff --git a/celltk/utils/tf_deepcell/tfutils.py b/celltk/utils/tf_deepcell/tfutils.py
deleted file mode 100644
index 6cceb62..0000000
--- a/celltk/utils/tf_deepcell/tfutils.py
+++ /dev/null
@@ -1,117 +0,0 @@
-
-import os
-import imp
-try:
-    from tensorflow.python.keras import backend
-    from tensorflow.python.keras.layers import Layer, Conv2D, MaxPooling2D
-    from tensorflow.python.keras.models import Sequential
-    from tensorflow.python.keras.models import load_model
-except:
-    from tensorflow.contrib.keras.python.keras.engine.topology import Layer
-    from tensorflow.contrib.keras.python.keras import backend
-    from tensorflow.contrib.keras.python.keras.layers import Conv2D, MaxPooling2D
-    from tensorflow.contrib.keras.python.keras.models import Sequential
-    from tensorflow.contrib.keras.python.keras.models import load_model
-from _dilated_pool import DilatedMaxPool2D
-import numpy as np
-from scipy.ndimage import imread as imread0
-import tifffile as tiff
-
-
-class Squeeze(Layer):
-    def __init__(self, output_dim=None, **kwargs):
-        self.output_dim = output_dim
-        super(Squeeze, self).__init__(**kwargs)
-
-    def call(self, x):
-        x = backend.squeeze(x, axis=2)
-        return backend.squeeze(x, axis=1)
-
-    def compute_output_shape(self, input_shape):
-        return (input_shape[0], input_shape[3])
-
-    # def get_config(self):
-    #     config = {'output_dim': self.output_dim}
-    #     base_config = super(Squeeze, self).get_config()
-    #     return dict(list(base_config.items()) + list(config.items()))
-
-
-
-def convert_model_patch2full(model):
-    """
-    """
-    dr = 1
-    new_model = Sequential()
-    for nl, layer in enumerate(model.layers):
-        if layer.name.startswith('squeeze_'):
-            continue
-        if isinstance(layer, MaxPooling2D):
-            newl = DilatedMaxPool2D(dilation_rate=dr)
-            newl = newl.from_config(layer.get_config())
-            newl.strides, newl.dilation_rate = (1, 1), dr
-            new_model.add(newl)
-            dr = dr * 2
-            continue
-        if isinstance(layer, Conv2D):
-            if not layer.kernel_size == (1, 1):
-                layer.dilation_rate = (dr, dr)
-                new_model.add(layer)
-            else:
-                newl = Conv2D(layer.filters, layer.kernel_size, input_shape=layer.input_shape[1:])
-                new_model.add(newl.from_config(layer.get_config()))
-        else:
-            new_model.add(layer)
-    return new_model
-
-
-def load_model_py(path):
-    if path.endswith('.py'):
-        fname = os.path.basename(path).split('.')[0]
-        module = imp.load_source(fname, path)
-        return module.model
-    elif path.endswith('.hdf5'):
-        return load_model(path, custom_objects={'Squeeze':Squeeze})
-
-
-
-def make_outputdir(output):
-    try:
-        os.makedirs(output)
-    except:
-        pass
-
-
-def imread_check_tiff(path):
-    img = imread0(path)
-    if img.dtype == 'object' or path.endswith('tif'):
-        img = tiff.imread(path)
-    return img
-
-
-def imread(path):
-    if isinstance(path, tuple) or isinstance(path, list):
-        st = []
-        for p in path:
-            st.append(imread_check_tiff(p))
-        img = np.dstack(st)
-        if img.shape[2] == 1:
-            np.squeeze(img, axis=2)
-        return img
-    else:
-        return imread_check_tiff(path)
-
-
-def parse_image_files(inputs):
-    if "/" not in inputs:
-        return (inputs, )
-    store = []
-    li = []
-    while inputs:
-        element = inputs.pop(0)
-        if element == "/":
-            store.append(li)
-            li = []
-        else:
-            li.append(element)
-    store.append(li)
-    return zip(*store)
\ No newline at end of file
diff --git a/celltk/utils/tf_deepcell/train.py b/celltk/utils/tf_deepcell/train.py
deleted file mode 100644
index f8dff1f..0000000
--- a/celltk/utils/tf_deepcell/train.py
+++ /dev/null
@@ -1,113 +0,0 @@
-from __future__ import division, print_function
-import os
-import numpy as np
-try:
-    from tensorflow.python.keras import optimizers, callbacks
-    from tensorflow.python.keras.preprocessing.image import ImageDataGenerator
-except:
-    from tensorflow.contrib.keras import optimizers, callbacks
-    from tensorflow.contrib.keras.python.keras.preprocessing.image import ImageDataGenerator
-from tfutils import imread
-from patches import extract_patches, pick_coords, pick_coords_list, extract_patch_list, _extract_patches, PatchDataGeneratorList
-from tfutils import load_model_py, make_outputdir
-from os.path import join
-from tfutils import parse_image_files
-
-FRAC_TEST = 0.1
-
-
-def define_callbacks(output, batch_size):
-    csv_logger = callbacks.CSVLogger(join(output, 'training.log'))
-    earlystop = callbacks.EarlyStopping(monitor='val_loss', patience=2)
-    tensorboard = callbacks.TensorBoard(batch_size=batch_size)
-    fpath = join(output, 'weights.{epoch:02d}-{loss:.2f}-{acc:.2f}-{val_loss:.2f}-{val_acc:.2f}.hdf5')
-    cp_cb = callbacks.ModelCheckpoint(filepath=fpath, monitor='val_loss', save_best_only=True)
-    return [csv_logger, earlystop, tensorboard, cp_cb]
-
-
-def train(image_list, labels_list, model_path, output, patchsize=61, nsamples=10000,
-          batch_size=32, nepochs=100, frac_test=FRAC_TEST):
-    assert np.bool(patchsize & 0x1)  # check if odd
-    model = load_model_py(model_path)
-    model.summary()
-
-    li_image, li_labels = [], []
-    for image_path, labels_path in zip(image_list, labels_list):
-        image, labels = imread(image_path), imread(labels_path).astype(np.uint8)
-        if image.ndim == 2:
-            image = np.expand_dims(image, -1)
-        elif image.ndim == 3:
-            image = np.moveaxis(image, 0, -1)
-        li_image.append(image)
-        li_labels.append(labels)
-
-    num_tests = int(nsamples * FRAC_TEST)
-    ecoords = pick_coords_list(nsamples, li_labels, patchsize, patchsize)
-    ecoords_tests, ecoords_train = ecoords[:num_tests], ecoords[num_tests:],
-    x_tests, y_tests = extract_patch_list(li_image, li_labels, ecoords_tests, patchsize, patchsize)
-    li_image = [np.expand_dims(i, 0) for i in li_image]
-
-    make_outputdir(output)
-    opt = optimizers.SGD(lr=0.01, decay=1e-6, momentum=0.9, nesterov=True)
-    model.compile(optimizer=opt, loss='sparse_categorical_crossentropy', metrics=['accuracy'])
-    callbacksets = define_callbacks(output, batch_size)
-
-    datagen = PatchDataGeneratorList(rotation_range=90, shear_range=0,
-                                     horizontal_flip=True, vertical_flip=True)
-    history = model.fit_generator(datagen.flow(li_image, li_labels, ecoords_train, patchsize, patchsize, batch_size=batch_size, shuffle=True),
-                                  steps_per_epoch=len(ecoords_train)/batch_size,
-                                  epochs=nepochs,
-                                  validation_data=(x_tests, y_tests),
-                                  validation_steps=len(ecoords_train)/batch_size,
-                                  callbacks=callbacksets)
-
-    score = model.evaluate(x_tests, y_tests, batch_size=batch_size)
-    print('score[loss, accuracy]:', score)
-    rec = dict(acc=history.history['acc'], val_acc=history.history['val_acc'],
-               loss=history.history['loss'], val_loss=history.history['val_loss'])
-    np.savez(join(output, 'records.npz'), **rec)
-
-    json_string = model.to_json()
-    open(join(output, 'cnn_model.json'), 'w').write(json_string)
-    model.save_weights(join(output, 'cnn_model_weights.hdf5'))
-    yaml_string = model.to_yaml()
-    open(join(output, 'cnn_model.yaml'), 'w').write(yaml_string)
-
-
-def _parse_command_line_args():
-    """
-    image:  Path to a tif or png file (e.g. data/nuc0.png).
-            To pass multiple image files (size can be varied), use syntax like
-            "-i im0.tif / im1.tif / im2.tif", and pass the same number of labels.
-    labels: (e.g. data/labels0.tif)
-    model:  path to a python file describing a model (e.g. data/tests_model.py)
-            or weights.*.hdf5 produced by ModelCheckPoint.
-            To resume training, pass the hdf5 file.
-    n:      A number of pixels for training. Use a large number (like 1,000,000)
-    batch:  Typically 128-512?
-    """
-
-    import argparse
-    parser = argparse.ArgumentParser(description='predict')
-    parser.add_argument('-i', '--image', help='image file path', nargs="*")
-    parser.add_argument('-l', '--labels', help='labels file path', nargs="*")
-    parser.add_argument('-m', '--model', help='path to python file or hdf5')
-    parser.add_argument('-o', '--output', default='.', help='output directory')
-    parser.add_argument('-n', '--nsamples', type=int, default=10000, help='number of samples')
-    parser.add_argument('-b', '--batch', type=int, default=64)
-    parser.add_argument('-e', '--epoch', type=int, default=50)
-    parser.add_argument('-p', '--patch', type=int, default=61,
-                        help='pixel size of image patches. make it odd')
-    return parser.parse_args()
-
-
-def _main():
-    args = _parse_command_line_args()
-    images = parse_image_files(args.image)[0]
-    labels = parse_image_files(args.labels)[0]
-    train(images, labels, args.model, args.output, args.patch,
-          args.nsamples, args.batch, args.epoch)
-
-
-if __name__ == "__main__":
-    _main()
diff --git a/celltk/utils/tfutils.py b/celltk/utils/tfutils.py
new file mode 100644
index 0000000..20b70b3
--- /dev/null
+++ b/celltk/utils/tfutils.py
@@ -0,0 +1,110 @@
+import os
+import imp
+import numpy as np
+from scipy.ndimage import imread as imread0
+import tifffile as tiff
+
+
+
+def conv_labels2dto3d(labels):
+    lbnums = np.unique(labels)
+    arr = np.zeros((labels.shape[0], labels.shape[1], len(lbnums)), np.uint8)
+    for i in lbnums:
+        arr[:, :, i] = labels == i
+    return arr
+
+
+def normalize(orig_img):
+    percentile = 99.9
+    high = np.percentile(orig_img, percentile)
+    low = np.percentile(orig_img, 100-percentile)
+    img = np.minimum(high, orig_img)
+    img = np.maximum(low, img)
+    img = (img - low) / (high - low)
+    return img
+
+
+def make_outputdir(output):
+    try:
+        os.makedirs(output)
+    except:
+        pass
+
+
+def imread_check_tiff(path):
+    img = imread0(path)
+    if img.dtype == 'object' or path.endswith('tif'):
+        img = tiff.imread(path)
+    return img
+
+
+def imread(path):
+    if isinstance(path, tuple) or isinstance(path, list):
+        st = []
+        for p in path:
+            st.append(imread_check_tiff(p))
+        img = np.dstack(st)
+        if img.shape[2] == 1:
+            np.squeeze(img, axis=2)
+        return img
+    else:
+        return imread_check_tiff(path)
+
+
+def parse_image_files(inputs):
+    if "/" not in inputs:
+        return (inputs, )
+    store = []
+    li = []
+    while inputs:
+        element = inputs.pop(0)
+        if element == "/":
+            store.append(li)
+            li = []
+        else:
+            li.append(element)
+    store.append(li)
+    return zip(*store)
+
+
+def pad_image(image):
+    # assumes the image is an np array of dimensions d0 x d1 x d2 x d3
+    # where d1 is height, d2 is width, and d3 is colors
+    # returns a list of padded image, hpadding (tuple), wpadding (tuple)
+
+    height = image.shape[1]
+    hdelta = 8 - (height % 8)
+    if (hdelta == 8):
+        hpadding = (0, 0)
+    elif (hdelta % 2) == 0:
+        hpadding = (int(hdelta/2.0), int(hdelta/2.0))
+    else:
+        hpadding = (int(hdelta/2.0), int(hdelta/2.0)+1)
+
+    width = image.shape[2]
+    wdelta = 8 - (width % 8)
+    if (wdelta == 8):
+        wpadding = (0, 0)
+    elif (wdelta % 2) == 0:
+        wpadding = (int(wdelta/2.0), int(wdelta/2.0))
+    else:
+        wpadding = (int(wdelta/2.0), int(wdelta/2.0)+1)
+    return [np.pad(image, ((0, 0), hpadding, wpadding, (0, 0)), 'constant', constant_values=0.0), hpadding, wpadding]
+
+
+def normalize_predictions(predictions):
+    # predictions is typically a 3 x h x w list
+    predictions = np.array(predictions)
+    num_rows = predictions.shape[1]
+    num_cols = predictions.shape[2]
+
+    for i in range(num_rows):
+        for j in range(num_cols):
+            prediction = predictions[:, i, j]
+            if (prediction.sum() != 0):
+                predictions[:, i, j] = prediction / prediction.sum()
+
+    return predictions
+
+
+
diff --git a/celltk/utils/tf_deepcell/predict.py b/celltk/utils/unet_predict.py
similarity index 52%
rename from celltk/utils/tf_deepcell/predict.py
rename to celltk/utils/unet_predict.py
index 9b3ede7..fa267cb 100644
--- a/celltk/utils/tf_deepcell/predict.py
+++ b/celltk/utils/unet_predict.py
@@ -3,40 +3,37 @@
 import numpy as np
 from os.path import join, basename, splitext
 from tfutils import imread
-try:
-    from tensorflow.python.keras import backend
-except:
-    from tensorflow.contrib.keras.python.keras import backend
-from tfutils import convert_model_patch2full, load_model_py, make_outputdir
+from tfutils import make_outputdir, normalize
+from tfutils import pad_image, normalize_predictions
 import tifffile as tiff
+import _model_builder
 
 
-def predict(img_path, model_path, weight_path):
+def predict(img_path, weight_path):
     x = imread(img_path)
+    x = normalize(x)
 
     if x.ndim == 2:
         x = np.expand_dims(x, -1)
     elif x.ndim == 3:
         x = np.moveaxis(x, 0, -1)
     x = np.expand_dims(x, 0)
+    num_colors = x.shape[-1]
 
-    model = load_model_py(model_path)
-    model = convert_model_patch2full(model)
-    model.load_weights(weight_path)
+    x, hpadding, wpadding = pad_image(x)
 
-    model.summary()
-    evaluate_model = backend.function(
-        [model.layers[0].input, backend.learning_phase()],
-        [model.layers[-1].output]
-        )
+    model = _model_builder.get_model(x.shape[1], x.shape[2], num_colors, activation=None)
+    model.load_weights(weight_path)
+    predictions = model.predict(x, batch_size=1)
+    predictions = [predictions[0, :, :, i] for i in range(predictions.shape[-1])]
 
-    cc = evaluate_model([x, 0])[0]
+    # resize predictions to match image dimensions (i.e. remove padding)
+    height = np.shape(predictions[0])[0]
+    width = np.shape(predictions[0])[1]
+    predictions = [p[hpadding[0]:height-hpadding[1], wpadding[0]:width-wpadding[1]] for p in predictions]
 
-    # from tensorflow.contrib.keras import optimizers
-    # opt = optimizers.SGD(lr=0.01, decay=1e-6, momentum=0.9, nesterov=True)
-    # model.compile(optimizer=opt, loss='sparse_categorical_crossentropy', metrics=['accuracy'])
-    # cc = model.predict(x)
-    return [cc[0, :, :, i] for i in range(cc.shape[-1])]
+    predictions = normalize_predictions(predictions)
+    return predictions
 
 
 def save_output(outputdir, images, pattern):
@@ -50,14 +47,13 @@ def _parse_command_line_args():
     parser = argparse.ArgumentParser(description='predict')
     parser.add_argument('-i', '--image', help='image file path')
     parser.add_argument('-w', '--weight', help='hdf5 file path')
-    parser.add_argument('-m', '--model', help='python file path with models')
     parser.add_argument('-o', '--output', default='.', help='output directory')
     return parser.parse_args()
 
 
 def _main():
     args = _parse_command_line_args()
-    images = predict(args.image, args.model, args.weight)
+    images = predict(args.image, args.weight)
     save_output(args.output, images, splitext(basename(args.image))[0])