diff --git a/calculate_cell_positions.py b/calculate_cell_positions.py
new file mode 100644
index 0000000..3ec14bf
--- /dev/null
+++ b/calculate_cell_positions.py
@@ -0,0 +1,67 @@
+import numpy as np
+import numpy.typing as npt
+
+TWO_FIVE_FIVE = 255
+
+
+def calculate_cell_positions(
+    image: npt.NDArray[np.float64], labelled_cells: npt.NDArray[np.uint16], type: bool
+) -> None:
+    """
+    This function calculates the min intensity position of each
+    labelled cell or the centroid position of each labelled region
+    """
+    no_cells = labelled_cells.max()
+
+    pos = np.zeros((no_cells, 2))
+
+    subtracted_image = TWO_FIVE_FIVE - image
+
+    for n in range(no_cells):
+
+        sy, sx = np.argwhere(labelled_cells == n).T
+
+        if type:
+            _place_at_lowest_int(pos, subtracted_image, sx, sy, n)
+
+        else:
+            _place_at_centroid(pos, image, sx, sy, n)
+
+            if ~np.isnan(pos[n, 1]) and labelled_cells[pos[n, 1], pos[n, 0]] != n:
+                # every so often the centroid is actually not in the label!
+                _place_at_lowest_int(pos, image, sx, sy, n)
+
+    return pos
+
+
+def _place_at_centroid(pos, image, sx, sy, n) -> None:
+    """
+    calculating the centroid from intensities
+    """
+    sum_x = 0
+    sum_y = 0
+    sum_I = 0
+
+    for m in range(len(sy)):
+        sum_x += sx[m] * image[sy[m], sx[m]]
+        sum_y += sy[m] * image[sy[m], sx[m]]
+        sum_I += image[sy[m], sx[m]]
+
+    pos[n, 1] = round(sum_y / sum_I)
+    pos[n, 0] = round(sum_x / sum_I)
+    return pos
+
+
+def _place_at_lowest_int(pos, image, sx, sy, n) -> None:
+    """
+    looking for the lowest intensity
+    """
+    val = TWO_FIVE_FIVE
+
+    for m in range(len(sy)):
+
+        if image[sy[m], sx[m]] < val:
+
+            val = image[sy[m], sx[m]]
+            pos[n, 2] = sy[m]
+            pos[n, 1] = sx[m]
diff --git a/example_segmentation.py b/example_segmentation.py
new file mode 100644
index 0000000..5b00dd3
--- /dev/null
+++ b/example_segmentation.py
@@ -0,0 +1,41 @@
+from pathlib import Path
+
+from scipy import io as sio
+from skimage.morphology import disk
+
+from segment_image import segment_image
+from utils import unlabel_poor_seeds_in_frame
+
+_file_location = Path(__file__).resolve()
+_matfile = _file_location.parent / "ProjIm.mat"
+
+
+def main() -> None:
+    """
+    implements the demo in the README in python
+    """
+    # load example image (Projected Drosophila Wing Disc - Ecad:GFP)
+    proj_im = sio.loadmat(_matfile)["ProjIm"]
+
+    # crop image for testing
+    cropped_image = proj_im[300:500, 400:600]
+
+    segmentation, seeds, labels = segment_image(
+        cropped_image,
+        sigma_1=0.5,
+        min_cell_size=2,
+        threshold=20,
+        merge_criteria=0.0,
+        sigma_3=0.5,
+        large_cell_size_thres=3000,
+        I_bound_max_pcnt=0.1,
+        show_feedback=True,
+    )
+
+
+if __name__ == "__main__":
+    se = disk(2)
+    Im = sio.loadmat(_file_location.parent / "Im.mat")["Im"]
+    CellSeeds = sio.loadmat(_file_location.parent / "CellSeeds.mat")["CellSeeds"]
+    CellLabels = sio.loadmat(_file_location.parent / "CellLabels.mat")["CellLabels"]
+    unlabel_poor_seeds_in_frame(Im, CellSeeds, CellLabels, se, 20, 0.5, 0.1),
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..81814bc
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,3 @@
+numpy==1.23.1
+scikit-image==0.19.3
+scipy==1.9.0
diff --git a/segment_image.py b/segment_image.py
new file mode 100644
index 0000000..50d33db
--- /dev/null
+++ b/segment_image.py
@@ -0,0 +1,86 @@
+import numpy as np
+import numpy.typing as npt
+from skimage.morphology import disk
+
+from utils import (
+    create_color_boundaries,
+    do_initial_seeding,
+    grow_cells_in_frame,
+    merge_seeds_from_labels,
+    neutralise_pts_not_under_label_in_frame,
+    unlabel_poor_seeds_in_frame,
+)
+
+
+def segment_image(
+    image: npt.NDArray[np.uint8],
+    *,
+    sigma_1: float,
+    min_cell_size: float,
+    threshold: int,
+    merge_criteria: float,
+    sigma_3: float,
+    large_cell_size_thres: float,
+    I_bound_max_pcnt: float,
+    show_feedback: bool = False
+) -> tuple[npt.NDArray[np.uint8], npt.NDArray[np.uint8], npt.NDArray[np.uint16]]:
+    """
+    Segments a single frame extracting the cell outlines.
+
+    Args:
+      image:
+        increasing intesity for membrane
+      sigma_1:
+        size px of gaussian for smoothing image [0+]
+      min_cell_size:
+        size px of smallest cell expected [0+]
+      threshold:
+        minimum value for membrane signal [0-255]
+      merge_criteria:
+        minimum ratio of low intensity pxs upon which to merge cells [0-1]
+      sigma_3:
+        size px of gaussian for smoothing image [0+]
+      large_cell_size_thres:
+        size px of largest cell expected [0+]
+      I_bound_max_pcnt:
+        minimum ratio for seed and membrane intensity [0-1]
+      show_feedback:
+        show feedback for segmentation
+
+    Returns:
+      segmented_image:
+        Im with seeds [255] and cell outlines [0]
+      cell_seeds:
+        Rescaled Image to fit the range [0,252]
+        253/254/255 are used for seed information
+      cell_labels:
+        bitmap of cells colored with 16bit id
+    """
+    # initialise
+    cell_seeds = np.zeros(image.shape, dtype=np.uint8)
+    # TODO: why using 16 bit for labels?
+    cell_labels = np.zeros(image.shape, dtype=np.uint16)
+
+    # TODO: check image casting
+    image = image.astype(float)
+    image *= 252 / image.max()
+    image = image.astype(np.uint8)
+
+    # structuring element, SE, used for morphological operations
+    se = disk(2)
+
+    do_initial_seeding(image, sigma_1, min_cell_size, threshold, large_cell_size_thres)
+
+    merge_seeds_from_labels(image, cell_seeds, cell_labels, se, merge_criteria, sigma_3)
+
+    grow_cells_in_frame(image, cell_seeds, sigma_3)
+
+    unlabel_poor_seeds_in_frame(
+        image, cell_seeds, cell_labels, se, threshold, sigma_3, I_bound_max_pcnt
+    )
+
+    neutralise_pts_not_under_label_in_frame(cell_seeds, cell_labels)
+
+    create_color_boundaries(image, cell_seeds, cell_labels)
+
+    return image, cell_seeds, cell_labels
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..8129364
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,8 @@
+[flake8]
+ignore = E203,E741,W503,W605
+max-complexity = 18
+max-line-length = 88
+select = B,B9,BLK,C,E,F,I,S,T4,W
+
+[mypy]
+plugins = numpy.typing.mypy_plugin
diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000..eef0ed5
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,392 @@
+import numpy as np
+import numpy.typing as npt
+from skimage.measure import regionprops
+from skimage.morphology import binary_dilation, binary_erosion
+
+from calculate_cell_positions import calculate_cell_positions
+
+SIGMA_BOUND = 0.01
+# TODO: work out what these are for and give better names
+TWO_FIVE_FIVE = 255
+TWO_FIVE_THREE = 253
+TWO_FIVE_TWO = 252
+
+
+def create_color_boundaries(
+    image: npt.NDArray[np.uint8],
+    cell_seeds: npt.NDArray[np.uint8],
+    cell_labels: npt.NDArray[np.uint16],
+) -> npt.NDArray[np.uint8]:
+    """
+    Generate final colored image (RGB) to represent the segmentation results
+    """
+    # given that every cell has a different label we can compute
+    # the boundaries by computing where the gradient changes
+    cell_labels = cell_labels.astype(float)
+    gy, gx = np.gradient(cell_labels)
+    cell_outlines = (cell_labels > 0) & ((gx**2 + gy**2) > 0)
+
+    cell_seeds_mask = cell_seeds > TWO_FIVE_THREE
+
+    image[cell_outlines] = 0
+    image[cell_seeds_mask] = TWO_FIVE_FIVE
+    return image
+
+
+def do_initial_seeding(
+    image: npt.NDArray[np.uint8],
+    sigma_1: float,
+    min_cell_size: float,
+    threshold: int,
+    large_cell_size_thres: float,
+) -> npt.NDArray[np.uint8]:
+    """
+    Find the initial cell seeds
+    """
+    # Create gaussian filter
+    if sigma_1 > SIGMA_BOUND:
+        f1 = _matlab_style_gauss2D(image.shape, sigma_1)
+
+        # Gaussian smoothing for the segmentation of individual cells
+        smoothed_image = np.fft.fftshift(
+            np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1))
+        ).real
+    else:
+        smoothed_image = image.astype(float)
+
+    smoothed_image /= smoothed_image.max() * TWO_FIVE_TWO
+
+    # Use external c-code to find initial seeds
+    initial_labelling = findcellsfromregiongrowing(
+        smoothed_image, min_cell_size, threshold
+    )
+
+    # set unallocated pixels to 0
+    initial_labelling[initial_labelling == 1] = 0
+
+    # Generate cell_labels from inital_labelling
+    cell_labels = initial_labelling.astype(np.uint16)
+
+    # eliminate very large areas
+    delabel_very_large_areas(cell_labels, large_cell_size_thres)
+
+    # Use true centre of cells as labels
+    centroids = np.round(calculate_cell_positions(smoothed_image, cell_labels, False))
+    centroids = centroids[~np.isnan(centroids.T[0])]
+    for n in range(len(centroids)):
+        smoothed_image[centroids[n, 1], centroids[n, 0]] = TWO_FIVE_FIVE
+
+    # cell_seeds contains the position of the true cell center.
+    return smoothed_image.astype(np.uint8)
+
+
+def grow_cells_in_frame(
+    image: npt.NDArray[np.uint8], cell_seeds: npt.NDArray[np.uint8], sigma_3: float
+) -> npt.NDArray[np.uint16]:
+    """
+    Growing cells from seeds TODO: add paramters in Name description!
+    """
+    # find labels
+    bw = (cell_seeds > TWO_FIVE_TWO).astype(float)
+
+    if sigma_3 > SIGMA_BOUND:
+        f1 = _matlab_style_gauss2D(image.shape, sigma_3)
+        smoothed_image = np.fft.fftshift(
+            np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1))
+        ).real
+    else:
+        smoothed_image = image.astype(float)
+
+    # mark labels on image
+    image_with_seeds = (smoothed_image).astype(float) * (1 - bw) + TWO_FIVE_FIVE * bw
+    return growcellsfromseeds3(image_with_seeds, TWO_FIVE_THREE).astype(np.uint16)
+
+
+def unlabel_poor_seeds_in_frame(
+    image: npt.NDArray[np.uint8],
+    cell_seeds: npt.NDArray[np.uint8],
+    cell_labels: npt.NDArray[np.uint16],
+    se: npt.NDArray[np.uint8],
+    threshold: int,
+    sigma_3: float,
+    I_bound_max_pcnt: float,
+) -> None:
+    """
+    Eliminate labels from seeds which have poor boundary intensity
+    """
+    L = cell_labels
+
+    if sigma_3 > SIGMA_BOUND:
+        f1 = _matlab_style_gauss2D(image.shape, sigma_3)
+        smoothed_image = np.fft.fftshift(
+            np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1))
+        ).real
+    else:
+        smoothed_image = image.astype(float)
+
+    # i.e. every cell is marked by one unique integer label
+    label_list = np.unique(L)
+    label_list = label_list[label_list != 0]
+    IBounds = np.zeros(len(label_list))
+    decisions = [0, 0, 0, 0, 0]
+
+    for c in range(len(label_list)):
+        mask = L == label_list[c]
+        cpy, cpx = np.argwhere(mask > 0).T
+        # find region of that label
+        minx = max(cpx.min() - 5, 0)
+        miny = max(cpy.min() - 5, 0)
+        maxx = min(cpx.max() + 5, image.shape[1] - 1)
+        maxy = min(cpy.max() + 5, image.shape[0] - 1)
+        # reduced to region of the boundary
+        reduced_mask = mask[miny : maxy + 1, minx : maxx + 1]
+        reduced_image = smoothed_image[miny : maxy + 1, minx : maxx + 1]
+        dilated_mask = binary_dilation(reduced_mask, se)
+        eroded_mask = binary_erosion(reduced_mask, se)
+        boundary_mask = dilated_mask ^ eroded_mask
+        boundary_intensities = reduced_image[boundary_mask > 0]
+        H = reduced_image[boundary_mask > 0]
+        I_bound = (boundary_intensities).mean()
+        IBounds[c] = I_bound
+
+        # cell seed information is retrieved as comparison
+        F2 = cell_seeds.copy()
+        F2[~mask] = 0
+        cpy, cpx = np.argwhere(F2 > TWO_FIVE_TWO).T
+        ICentre = smoothed_image[cpy, cpx][0]
+
+        I_bound_max = TWO_FIVE_FIVE * I_bound_max_pcnt
+
+        # Figure out which conditions make the label invalid
+        # 1. I_bound_max, gives the Lower bound to the mean intensity
+        # 1.b condition upon that the cell seed has less
+        #     than 20% intensity difference to the mean
+        #   => If the cell boundary is low and not very different
+        #      from the seed, cancel the region
+        first_condition = (I_bound < I_bound_max) & (I_bound / ICentre < 1.2)
+        # 2. W/o (1.b) the lower bound is reduced by ~17% (1 - 0.833) to be decisive
+        second_condition = I_bound < I_bound_max * 5 / 6
+        # 3. If the minimum retrieved in the boundary mask is 0 (dangerous!)
+        third_condition = np.min(boundary_intensities) == 0
+        # 4. If the amount of low intensity signal (i.e. < 20) is more than 10%
+        fourth_condition = np.sum(H < threshold) / len(H) > 0.1
+        if first_condition | second_condition | third_condition | fourth_condition:
+
+            # The label is cancelled (inverted mask multiplication.)
+            cell_labels *= (mask == 0).astype(np.uint16)
+
+            # record the removal decisions
+            if first_condition:
+                decisions[1] += 1
+            elif second_condition:
+                decisions[2] += 1
+            elif third_condition:
+                decisions[3] += 1
+            elif fourth_condition:
+                decisions[4] += 1
+            else:
+                # should not happen
+                decisions[5] += 1
+
+
+def delabel_very_large_areas(
+    cell_labels: npt.NDArray[np.uint16], large_cell_size_thres: float
+) -> npt.NDArray[np.uint16]:
+    """
+    remove cells which are bigger than large_cell_size_thres
+    """
+    As = [r.area for r in regionprops(cell_labels, cache=False)]
+    ls = np.unique(cell_labels)
+    for l in range(len(ls)):
+        if l == 0:
+            continue
+
+        A = As[l]
+        if A > large_cell_size_thres:
+            cell_labels[cell_labels == l] = 0
+
+    return cell_labels
+
+
+def merge_seeds_from_labels(
+    image: npt.NDArray[np.uint8],
+    cell_seeds: npt.NDArray[np.uint8],
+    cell_labels: npt.NDArray[np.uint16],
+    se: npt.NDArray[np.uint8],
+    merge_criteria: float,
+    sigma_3: float,
+) -> None:
+    """
+    Remove initial cell regions which touch & whose boundary is insufficient
+    """
+    # smoothing
+    if sigma_3 > SIGMA_BOUND:
+        f1 = _matlab_style_gauss2D(image.shape, sigma_3)
+        smoothed_image = np.fft.fftshift(
+            np.fft.ifft2(np.fft.fft2(image) * np.fft.fft2(f1))
+        ).real
+    else:
+        smoothed_image = image.astype(float)
+
+    label_list = np.unique(cell_labels)
+    label_list = label_list[label_list != 0]
+    c = 1
+
+    merge_intensity_distro: list = []
+    merge_decisions = 0
+
+    # loop over labels
+    while True:
+        label_mask = cell_labels == label_list[c]
+        label = label_list[c]
+
+        cpy, cpx = np.argwhere(label_mask > 0).T
+
+        # find region of that label
+        minx = np.min(cpx)
+        maxx = np.max(cpx)
+        miny = np.min(cpy)
+        maxy = np.max(cpy)
+        minx = np.max(minx - 5, 1)
+        miny = np.max(miny - 5, 1)
+        maxx = np.min(maxx + 5, image.shape[1])
+        maxy = np.min(maxy + 5, image.shape[0])
+
+        # reduce data to that region
+        reduced_label_mask = label_mask[miny:maxy, minx:maxx]
+        reduced_image = smoothed_image[miny:maxy, minx:maxx]
+        reduced_labels = cell_labels[miny:maxy, minx:maxx]
+
+        # now find boundaries
+        dilated_mask = binary_dilation(reduced_label_mask, se)
+        eroded_mask = binary_erosion(reduced_label_mask, se)
+        border_mask = dilated_mask - eroded_mask
+        border_intensities = reduced_image[border_mask > 0]
+        central_intensity = reduced_image[eroded_mask > 0]
+
+        F2 = cell_seeds
+        F2[~label_mask] = 0
+        cpy, cpx = np.argwhere(F2 > TWO_FIVE_THREE).T
+        ICentre = smoothed_image[cpy, cpx]
+
+        background_std = np.std(central_intensity.astype(float))
+
+        # get labels of surrounding cells (neighbours)
+        neighbour_labels = np.unique(reduced_labels[dilated_mask > 0])
+        neighbour_labels = neighbour_labels[neighbour_labels != label]
+
+        low_intensity_ratios: list = []
+        for i in range(len(neighbour_labels)):
+            neighb_label = neighbour_labels(i)
+            neighbor_border = dilated_mask
+            # slice of neighbour around cell
+            neighbor_border[reduced_labels != neighb_label] = 0
+            cell_border = binary_dilation(neighbor_border, se)
+            # slice of cell closest to neighbour
+            cell_border[reduced_labels != label] = 0
+
+            # combination of both creating boundary region
+            joint_border = (cell_border + neighbor_border) > 0
+            border_intensities = reduced_image
+            # intensities at boundary
+            border_intensities[~joint_border] = 0
+
+            # average number of points in boundary where intensity is
+            # of low quality (dodgy)
+            low_intensity_threshold = ICentre + (background_std / 2)
+            low_intensity_pixels = (
+                border_intensities[joint_border] < low_intensity_threshold
+            )
+
+            low_intensity_ratio = np.sum(low_intensity_pixels) / np.shape(
+                border_intensities[joint_border]
+            )
+
+            low_intensity_ratios = [low_intensity_ratios, low_intensity_ratio]
+
+        # Find out which is border with the lowest intensity ratio
+        worst_intensity_ratio, worst_neighbor_index = np.max(low_intensity_ratios)
+        neighb_label = neighbour_labels(worst_neighbor_index)
+
+        # if the label value is of poor quality, then recursively check the merge
+        # criteria in order to add it as a potential label in the label set.
+
+        merge_intensity_distro = [merge_intensity_distro, worst_intensity_ratio]
+
+        if (worst_intensity_ratio > merge_criteria) & label != 0 & neighb_label != 0:
+
+            merge_labels(cell_seeds, cell_labels, label, neighb_label)
+            label_list = np.unique(cell_labels)
+            label_list = label_list[label_list != 0]
+            c -= 1
+            # reanalyze the same cell for more possible mergings
+            merge_decisions += 1
+
+        c += 1
+
+        # Condition to break the while cycle -> as soon as all the
+        # labels are processed, then exit
+        if c > len(label_list):
+            break
+
+
+def merge_labels(
+    cell_seeds: npt.NDArray[np.uint8], cell_labels: npt.NDArray[np.uint16], l1, l2
+) -> None:
+    """ """
+    Cl = cell_labels
+    Il = cell_seeds
+    m1 = Cl == l1
+    m2 = Cl == l2
+    Il1 = Il
+    Il1[~m1] = 0
+    Il2 = Il
+    Il2[~m2] = 0
+    cpy1, cpx1 = np.argwhere(Il1 > TWO_FIVE_THREE).T
+    cpy2, cpx2 = np.argwhere(Il2 > TWO_FIVE_THREE).T
+    cpx = round((cpx1 + cpx2) / 2)
+    cpy = round((cpy1 + cpy2) / 2)
+
+    # background level
+    cell_seeds[cpy1, cpx1] = 20
+    cell_seeds[cpy2, cpx2] = 20
+    if (cell_labels[cpy, cpx] == l1) | (cell_labels[cpy, cpx] == l2):
+        cell_seeds[cpy, cpx] = TWO_FIVE_FIVE
+    elif np.sum(m1) > np.sum(m2):
+        cell_seeds[cpy1, cpx1] = TWO_FIVE_FIVE
+    else:
+        cell_seeds[cpy2, cpx2] = TWO_FIVE_FIVE
+
+    Cl[m2] = l1
+    cell_labels = Cl
+
+
+def neutralise_pts_not_under_label_in_frame(
+    cell_seeds: npt.NDArray[np.uint8],
+    cell_labels: npt.NDArray[np.uint16],
+) -> npt.NDArray[np.uint8]:
+    """
+    Seeds whose label has been eliminated are converted to neutral seeds
+
+    the idea here is to set seeds not labelled to a value
+    ie invisible to retracking (and to growing, caution!)
+    """
+    cell_seeds_copy = cell_seeds.copy()
+    cell_seeds_copy[cell_labels != 0] = 0
+    cell_seeds[cell_seeds_copy > TWO_FIVE_TWO] = TWO_FIVE_THREE
+    return cell_seeds
+
+
+def _matlab_style_gauss2D(shape: tuple[int, ...], sigma: float):
+    """
+    2D gaussian mask - should give the same result as MATLAB's
+    fspecial('gaussian', [shape], [sigma])
+    """
+    m, n = [(ss - 1) / 2 for ss in shape]
+    y, x = np.ogrid[-m : m + 1, -n : n + 1]
+    h = np.exp(-(x * x + y * y) / (2 * sigma * sigma))
+    h[h < np.finfo(h.dtype).eps * h.max()] = 0
+    sumh = h.sum()
+    if sumh != 0:
+        h /= sumh
+    return h