diff --git a/.github/workflows/docker-publish.yml b/.github/workflows/docker-publish.yml
index 8a63bef..6b0342b 100644
--- a/.github/workflows/docker-publish.yml
+++ b/.github/workflows/docker-publish.yml
@@ -35,11 +35,11 @@ jobs:
 
       # Install the cosign tool except on PR
       # https://github.com/sigstore/cosign-installer
-      - name: Install cosign
+      - name: Install Cosign
         if: github.event_name != 'pull_request'
-        uses: sigstore/cosign-installer@f3c664df7af409cb4873aa5068053ba9d61a57b6 #v2.6.0
+        uses: sigstore/cosign-installer@v3.1.1
         with:
-          cosign-release: 'v1.11.0'
+          cosign-release: 'v2.2.1'
           
       # Workaround: https://github.com/docker/build-push-action/issues/461
       - name: Setup Docker buildx
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..5b25f54
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,23 @@
+# schapirolabor/background_subtraction: Changelog
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## v0.4.1 - [2023.11.21]
+
+Complete rework of Backsub to include Palom's pyramid writer (https://github.com/labsyspharm/palom).
+Added dask array chunking and delayed execution for subtraction that happenes while the output pyramidal `ome.tif` is being created.
+Added `CHANGELOG.md`.
+
+### `Added`
+- `--chunk-size` parameter
+- Palom's pyramid writer
+
+### `Fixed`
+- Fixed issue with RAM inefficiency - reworked Backsub.
+
+### `Removed`
+- `--pyramid` tag introduced in v0.3.4, for smaller images, a smaller tile size should be specified now.
+
+
+I did not keep a changelog before version v0.4.1.
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
index 5cb0955..91574fc 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -1,7 +1,13 @@
-FROM mambaorg/micromamba:0.26.0
-COPY --chown=$MAMBA_USER:$MAMBA_USER environment.yml /tmp/env.yaml
-RUN micromamba install -y -n base -f /tmp/env.yaml && \
-    micromamba clean --all --yes
-ENV PATH="${PATH}:/opt/conda/bin"
+FROM continuumio/miniconda3
+
+COPY environment.yml .
+RUN apt-get update -qq && apt-get install -y \
+    build-essential \
+    ffmpeg \
+    libsm6 \
+    libxext6
+
+RUN conda env create -f environment.yml
+ENV PATH="/opt/conda/envs/backsub/bin:$PATH"
 WORKDIR /background_subtraction
-COPY . .
+COPY . .
\ No newline at end of file
diff --git a/README.md b/README.md
index 2428a95..03b7f18 100644
--- a/README.md
+++ b/README.md
@@ -13,13 +13,13 @@ Marker<sub>*corrected*</sub> = Marker<sub>*raw*</sub> - Background / Exposure<su
 
 ## Usage 
 
-### Markers file
-
-#### Versions v0.3.0 and newer:
 The `markers.csv` file which gives details about the channels needs to contain the following columns: "marker_name", "background" and "exposure". An exemplary [markers.csv](https://github.com/SchapiroLabor/Background_subtraction/blob/main/example/markers.csv) file is given. The "marker_name" column should indicate the marker for the acquired channel and all values should be unique. The "background" column should indicate the marker name of the channel which needs to be subtracted. This value must match the "marker_name" value of the background channel. The "exposure" column should contain the exposure time used for channel acquisition, and the measure unit should be consistent across the column. Exposure time is used for scaling the value of the background to be comparable to the processed channel. The "remove" column should contain logical `TRUE` values for channels which should be exluded in the output image.
 
+### Versions v0.4.1 and newer:
+The script has been rewritten to perform channel subtraction in a RAM-efficient manner - updating is highly recommended. If the output file is much bigger than expected, adjust the `--tile-size` parameter to a smaller value (e.g `512`). Changing the `--chunk-size` parameter may affect performance (lower values increase execution time, higher values increase RAM usage).
+
 
-#### Versions v0.2.0 and older:
+### Versions v0.2.0 and older:
 The `markers.csv` file which gives details about the channels needs to contain the following columns: "Filter", "background" and "exposure". An exemplary [markers_old.csv](https://github.com/SchapiroLabor/Background_subtraction/files/9549686/markers.csv) file is given. The "Filter" column should specify the Filter used when acquiring images. If different stains are aquired with the same filter, the *exact same value* needs to be written (including background) as it is used for determining which background channel should be subtracted. The "background" column should contain logical `TRUE` values for channels which represent autofluorescence. The "exposure" column should contain the exposure time used for channel acquisition, and the measure unit should be consistent across the column. Exposure time is used for scaling the value of the background to be comparable to the processed channel.
 
 
@@ -33,8 +33,8 @@ The script requires four inputs:
 Optional inputs:
 * `--pixel-size` to specify the pixel size of the input image (default: `1.0`), if not specified, the pixel size will be read from the metadata of the input image.
 * `--version` to print version and exit (added in v0.3.4)
-* `--pyramid` to create a pyramidal output image (default: `True`) (added in v0.3.4)
 * `--tile-size` to specify the tile size for the pyramidal output image (default: `1024`) (added in v0.3.4)
+* `--chunk-size` to specify chunk size for delayed calculation execution (default: `5000`) - lower values result in higher execution time, higher values in higher RAM usage
 
 
 ### Output
diff --git a/background_sub.py b/background_sub.py
index 5679c4d..31a5b93 100644
--- a/background_sub.py
+++ b/background_sub.py
@@ -1,30 +1,353 @@
-from __future__ import print_function, division
-from multiprocessing.spawn import import_main_path
-import sys
-import copy
 import argparse
+import pathlib
+
+import ome_types
+import palom.pyramid
+import palom.reader
+import pandas as pd
 import numpy as np
+import copy
+import dask.array as da
+import numexpr as ne
+import math
+import cv2
 import tifffile
+import tqdm
 import zarr
-import skimage.transform
-from aicsimageio import aics_image as AI
-import pandas as pd
-import numexpr as ne
-from ome_types import from_tiff, to_xml
-from os.path import abspath
+from loguru import logger
 from argparse import ArgumentParser as AP
+from os.path import abspath
 import time
-import dask
-# This API is apparently changing in skimage 1.0 but it's not clear to
-# me what the replacement will be, if any. We'll explicitly import
-# this so it will break loudly if someone tries this with skimage 1.0.
-try:
-    from skimage.util.dtype import _convert as dtype_convert
-except ImportError:
-    from skimage.util.dtype import convert as dtype_convert
+### Using the palom pyramidal writer by Yu-An Chen: https://github.com/labsyspharm/palom/blob/main/palom/pyramid.py
+### Adapted for this use case - only one image can be processed:
+
+class PyramidSetting:
+
+    def __init__(
+        self,
+        downscale_factor=2,
+        tile_size=1024,
+        max_pyramid_img_size=1024
+    ):
+        self.downscale_factor = downscale_factor
+        self.tile_size = tile_size
+        self.max_pyramid_img_size = max_pyramid_img_size
+
+    def tile_shapes(self, base_shape):
+        shapes = np.array(self.pyramid_shapes(base_shape))
+        n_rows_n_cols = np.ceil(shapes / self.tile_size)
+        tile_shapes = np.ceil(shapes / n_rows_n_cols / 16) * 16
+        return [tuple(map(int, s)) for s in tile_shapes]
+
+    def pyramid_shapes(self, base_shape):
+        num_levels = self.num_levels(base_shape)
+        factors = self.downscale_factor ** np.arange(num_levels)
+        shapes = np.ceil(np.array(base_shape) / factors[:,None])
+        return [tuple(map(int, s)) for s in shapes]
+
+    def num_levels(self, base_shape):
+        factor = max(base_shape) / self.max_pyramid_img_size
+        
+        return math.ceil(math.log(max(factor, 1), self.downscale_factor))+1
+
+def format_channel_names(num_channels_each_mosaic, channel_names):
+    '''
+    format_channel_names(
+        [1, 2, 3, 4, 5], ['x', 'x', ['x'], ['x', 'x']]
+    )
+    >>> [
+        'x_1',
+        'x_2',
+        'x_3',
+        'x_4',
+        'x_5',
+        'x_6',
+        'Mosaic 4_1',
+        'Mosaic 4_2',
+        'Mosaic 4_3',
+        'Mosaic 4_4'
+    ]
+    '''
+    matched_channel_names = []
+    for idx, (n, c) in enumerate( zip(channel_names, num_channels_each_mosaic*num_channels_each_mosaic[0]) ):
+        c=1
+        nl = n
+        if type(n) == str:
+            nl = [n]
+        if len(nl) == 1:
+            nl = nl
+        matched_channel_names.append(nl)
+    return make_unique_str(
+        [n for l in matched_channel_names for n in l]
+    )
+
+def make_unique_str(str_list):
+    if len(set(str_list)) == len(str_list):
+        return str_list
+    else:
+        max_length = max([len(s) for s in str_list])
+        str_np = np.array(str_list, dtype=np.dtype(('U', max_length+10)))
+        unique, counts = np.unique(str_np, return_counts=True)
+        has_duplicate = unique[counts > 1]
+        for n in has_duplicate:
+            suffixes = [
+                f"_{i}"
+                for i in range(1, (str_np == n).sum()+1)
+            ]
+            str_np[str_np == n] = np.char.add(n, suffixes)
+    return make_unique_str(list(str_np))
+
+def normalize_mosaics(mosaics):
+    dtypes = set(m.dtype for m in mosaics)
+    if any([np.issubdtype(d, np.floating) for d in dtypes]):
+        max_dtype = np.dtype(np.float32)
+    else:
+        max_dtype = max(dtypes)
+    normalized = []
+    for m in mosaics:
+        assert m.ndim == 2 or m.ndim == 3
+        if m.ndim == 2:
+            m = m[np.newaxis, :]
+        normalized.append(m.astype(max_dtype, copy=False))
+    return normalized
+
+
+def write_pyramid(
+    mosaics,
+    output_path,
+    pixel_size=1,
+    channel_names=None,
+    verbose=True,
+    downscale_factor=4,
+    compression=None,
+    is_mask=False,
+    tile_size=None,
+    save_RAM=False,
+    kwargs_tifffile=None
+):
+    mosaics = normalize_mosaics(mosaics)
+    ref_m = mosaics[0]
+    path = output_path
+    num_channels = count_num_channels(mosaics)
+    base_shape = ref_m.shape[1:3]
+    assert int(downscale_factor) == downscale_factor
+    assert downscale_factor < min(base_shape)
+    pyramid_setting = PyramidSetting(
+        downscale_factor=int(downscale_factor),
+        tile_size=max(ref_m.chunksize)
+    )
+    num_levels = pyramid_setting.num_levels(base_shape)
+    tile_shapes = pyramid_setting.tile_shapes(base_shape)
+    shapes = pyramid_setting.pyramid_shapes(base_shape)
+
+    if tile_size is not None:
+        assert tile_size % 16 == 0, (
+            f"tile_size must be None or multiples of 16, not {tile_size}"
+        )
+        tile_shapes = [(tile_size, tile_size)] * num_levels
+    print(tile_shapes)
+    dtype = ref_m.dtype
+
+    software = f'backsub {_version}'
+    metadata = {
+        'Creator': software,
+        'Pixels': {
+            'PhysicalSizeX': pixel_size, 'PhysicalSizeXUnit': '\u00b5m',
+            'PhysicalSizeY': pixel_size, 'PhysicalSizeYUnit': '\u00b5m'
+        },
+    }
+
+    if channel_names is not None:
+        num_channels_each_mosaic = [
+            count_num_channels([m])
+            for m in mosaics
+        ]
+        names = format_channel_names(num_channels_each_mosaic, channel_names)
+        if len(names) == num_channels:
+            metadata.update({
+                'Channel': {'Name': names},
+            })
+
+    logger.info(f"Writing to {path}")
+    with tifffile.TiffWriter(path, bigtiff=True) as tif:
+        kwargs = dict(
+            metadata=metadata,
+            software=software,
+            compression=compression
+        )
+        if kwargs_tifffile is None:
+            kwargs_tifffile = {}
+        tif.write(
+            data=tile_from_combined_mosaics(
+                mosaics, tile_shape=tile_shapes[0], save_RAM=save_RAM
+            ),
+            shape=(num_channels, *shapes[0]),
+            subifds=int(num_levels - 1),
+            dtype=dtype,
+            tile=tile_shapes[0],
+            **{**kwargs, **kwargs_tifffile}
+        )
+        logger.info('Generating pyramid')
+        for level, (shape, tile_shape) in enumerate(
+            zip(shapes[1:], tile_shapes[1:])
+        ):
+            if verbose:
+                logger.info(f"    Level {level+1} ({shape[0]} x {shape[1]})")
+            tif.write(
+                data=tile_from_pyramid(
+                    path,
+                    num_channels,
+                    tile_shape=tile_shape,
+                    downscale_factor=downscale_factor,
+                    level=level,
+                    is_mask=is_mask,
+                    save_RAM=save_RAM
+                ),
+                shape=(num_channels, *shape),
+                subfiletype=1,
+                dtype=dtype,
+                tile=tile_shape,
+                **{
+                    **dict(compression=compression),
+                    **kwargs_tifffile
+                }
+            )
+
 
+def count_num_channels(imgs):
+    for img in imgs:
+        assert img.ndim == 2 or img.ndim == 3
+    return sum([
+        1 if img.ndim == 2 else img.shape[0]
+        for img in imgs
+    ])
+   
+
+def tile_from_combined_mosaics(mosaics, tile_shape, save_RAM=False):
+    num_rows, num_cols = mosaics[0].shape[1:3]
+    h, w = tile_shape
+    n = len(mosaics)
+    for idx, m in enumerate(mosaics):
+        for cidx, c in enumerate(m):
+            # the performance is heavily degraded without pre-computing the
+            # mosaic channel
+            with tqdm.dask.TqdmCallback(
+                ascii=True,
+                desc=(
+                    f"Assembling mosaic {idx+1:2}/{n:2} (channel"
+                    f" {cidx+1:2}/{m.shape[0]:2})"
+                ),
+            ):
+                c = da_to_zarr(c) if save_RAM else c.compute()
+            for y in range(0, num_rows, h):
+                for x in range(0, num_cols, w):
+                    yield np.array(c[y:y+h, x:x+w])
+                    # yield m[y:y+h, x:x+w].copy().compute()
+            c = None
+
+def tile_from_pyramid(
+    path,
+    num_channels,
+    tile_shape,
+    downscale_factor=2,
+    level=0,
+    is_mask=False,
+    save_RAM=False
+):
+    # workaround progress bar
+    # https://forum.image.sc/t/tifffile-ome-tiff-generation-is-taking-too-much-ram/41865/26
+    pbar = tqdm.tqdm(total=num_channels, ascii=True, desc=f'Processing channel')
+    for c in range(num_channels):
+        img = da.from_zarr(zarr.open(tifffile.imread(
+            path, series=0, level=level, aszarr=True
+        )))
+        if img.ndim == 2:
+            img = img.reshape(1, *img.shape)
+        img = img[c]
+        # read using key seems to generate a RAM spike
+        # img = tifffile.imread(path, series=0, level=level, key=c)
+        if not is_mask:
+            img = img.map_blocks(
+                cv2.blur,
+                ksize=(downscale_factor, downscale_factor), anchor=(0, 0)
+            )
+        img = da_to_zarr(img) if save_RAM else img.compute()
+        num_rows, num_columns = img.shape
+        h, w = tile_shape
+        h *= downscale_factor
+        w *= downscale_factor
+        last_c = range(num_channels)[-1]
+        last_y = range(0, num_rows, h)[-1]
+        last_x = range(0, num_columns, w)[-1]
+        for y in range(0, num_rows, h):
+            for x in range(0, num_columns, w):
+                if (y == last_y) & (x == last_x):
+                    pbar.update(1)
+                    if c == last_c:
+                        pbar.close()
+                yield np.array(img[y:y+h:downscale_factor, x:x+w:downscale_factor])
+        # setting img to None seems necessary to prevent RAM spike
+        img = None
+
+
+def da_to_zarr(da_img, zarr_store=None, num_workers=None, out_shape=None, chunks=None):
+    if zarr_store is None:
+        if out_shape is None:
+            out_shape = da_img.shape
+        if chunks is None:
+            chunks = da_img.chunksize
+        zarr_store = zarr.create(
+            out_shape,
+            chunks=chunks,
+            dtype=da_img.dtype,
+            overwrite=True
+        )
+    da_img.to_zarr(zarr_store, compute=False).compute(
+        num_workers=num_workers
+    )
+    return zarr_store
+
+
+def process_markers(markers):
+    markers['ind'] = range(0, len(markers))
+    if 'remove' not in markers:
+        markers['remove'] = ["False" for i in range(len(markers))]
+    else:
+        markers['remove'] = markers['remove'] == True
+
+    markers['keep'] = markers['remove'] == False
+
+    markers = markers.drop(columns=['remove'])
+
+    return markers
+
+def detect_pixel_size(img_path,pixel_size=None):
+    if pixel_size is None:
+        print('Pixel size overwrite not specified')
+        try:
+            metadata = ome_types.from_tiff(img_path)
+            pixel_size = metadata.images[0].pixels.physical_size_x
+        except Exception as err:
+            print(err)
+            print('Pixel size detection using ome-types failed')
+            pixel_size = None
+    return pixel_size
+
+def scale_background(background_channel, scalar):
+    background_channel = np.rint(ne.evaluate("background_channel * scalar"))
+    return np.where(background_channel>65535,65535,background_channel.astype(np.uint16))
+
+def subtract(channel_to_process, background):
+    return np.where(channel_to_process<background, 0, channel_to_process-background)
+
+def subtract_channels(channel_to_process, background_channel, scalar, chunk_size):
+    channel_to_process = channel_to_process.rechunk(chunk_size)
+    background_channel = background_channel.rechunk(chunk_size)
+    background = da.map_blocks(scale_background, background_channel, scalar)
+    output = da.map_blocks(subtract, channel_to_process, background)
+    background = None
+    return output
 
-# arg parser
 def get_args():
     # Script description
     description="""Subtracts background - Lunaphore platform"""
@@ -36,10 +359,11 @@ def get_args():
     inputs = parser.add_argument_group(title="Required Input", description="Path to required input file")
     inputs.add_argument("-r", "--root", dest="root", action="store", required=True, help="File path to root image file.")
     inputs.add_argument("-m", "--markers", dest="markers", action="store", required=True, help="File path to required markers.csv file")
-    inputs.add_argument("--pixel-size", metavar="SIZE", dest = "pixel_size", type=float, default = None, action = "store",help="pixel size in microns; default is 1.0")
-    inputs.add_argument("--pyramid", dest="pyramid", required=False, default=True, help="Should output be pyramidal?")
-    inputs.add_argument("--tile-size", dest="tile_size", required=False, default=1024, help="Tile size for pyramid generation")
-    inputs.add_argument("--version", action="version", version="v0.3.4")
+    inputs.add_argument("--pixel-size", metavar="SIZE", dest = "pixel_size", type=float, default = None, action = "store",help="Pixel size in microns")
+    inputs.add_argument("--tile-size", dest="tile_size", required=False, type=int, default=1024, help="Tile size for pyramid generation")
+    #inputs.add_argument("--save-ram", dest="ram", required=False, default=False, help="Save RAM during pyramid generation")
+    inputs.add_argument("--chunk-size", dest="chunksize", required=False, type=int, default=5000, help="Chunk size for dask array (e.g for chunksize 1000, the image will be split into 1000x1000 chunks)")
+    inputs.add_argument("--version", action="version", version="v0.4.1")
     
     outputs = parser.add_argument_group(title="Output", description="Path to output file")
     outputs.add_argument("-o", "--output", dest="output", action="store", required=True, help="Path to output file")
@@ -55,253 +379,49 @@ def get_args():
 
     return arg
 
-def preduce(coords, img_in, img_out):
-    print(img_in.dtype)
-    (iy1, ix1), (iy2, ix2) = coords
-    (oy1, ox1), (oy2, ox2) = np.array(coords) // 2
-    tile = skimage.img_as_float32(img_in[iy1:iy2, ix1:ix2])
-    tile = skimage.transform.downscale_local_mean(tile, (2, 2))
-    tile = dtype_convert(tile, 'uint16')
-    #tile = dtype_convert(tile, img_in.dtype)
-    img_out[oy1:oy2, ox1:ox2] = tile
+def main(args):
+    _version = 'v0.4.1'
+    in_path = args.root
 
-def format_shape(shape):
-    return "%dx%d" % (shape[1], shape[0])
+    out_path = pathlib.Path(args.output)
+    # pixel data is read into RAM lazily, cannot overwrite input file
+    assert out_path != in_path
 
-def process_markers(markers):
-    # add column with starting indices (which match the image channel indices)
-    # this should be removed soon
-    markers['ind'] = range(0, len(markers))
+    # Detect pixel size in ome-xml
+    pixel_size = detect_pixel_size(in_path, args.pixel_size)
+    if pixel_size is None: pixel_size = 1
 
-    # if the 'remove' column is not specified, all channels are kept. If it is 
-    # present, it is converted to a boolean indicating which channels should be removed
-    if 'remove' not in markers:
-        markers['remove'] = ["False" for i in range(len(markers))]
-    else:
-        markers['remove'] = markers['remove'] == True
-    # invert the markers.remove column to indicate which columns to keep
-    markers['remove'] = markers['remove'] == False
+    markers = process_markers(pd.read_csv(args.markers))
 
-    return markers
+    # Use palom to lazily read the input image
+    readers = palom.reader.OmePyramidReader(in_path)
+    mosaic = readers.pyramid[0]
+    mosaic_out = copy.copy(mosaic)
+    chunk_size = (args.chunksize,args.chunksize)
+    mosaic_out = mosaic_out.rechunk((1,args.chunksize,args.chunksize))
 
-def process_metadata(metadata, markers):
-    try:
-        metadata.screens[0].reagents = [metadata.screens[0].reagents[i] for i in markers[markers.remove == True].ind]
-    except:
-        pass
-    try:
-        metadata.structured_annotations = [metadata.structured_annotations[i] for i in markers[markers.remove == True].ind]
-    except:
-        pass
-    # these two are required
-    metadata.images[0].pixels.size_c = len(markers[markers.remove == True])
-    metadata.images[0].pixels.channels = [metadata.images[0].pixels.channels[i] for i in markers[markers.remove == True].ind]
-    try:
-        metadata.images[0].pixels.planes = [metadata.images[0].pixels.planes[i] for i in markers[markers.remove == True].ind]
-    except:
-        pass
-    try:
-        metadata.images[0].pixels.tiff_data_blocks[0].plane_count = sum(markers.remove == True)
-    except:
-        pass
-    return metadata
-
-
-# NaN values return True for the statement below in this version of Python. Did not use math.isnan() since the values
-# are strings if present
-def isNaN(x):
-    return x != x
-
-def subtract_channel(image, markers, channel, background_marker, output):
-    scalar = markers[markers.ind == channel].exposure.values / background_marker.exposure.values
-    
-    # create temporary dataframe which will store the multiplied background rounded up to nearest integer
-    # [0] at the end needed to get [x, y] shape, and not have [1, x, y]
-    back = copy.copy(image[background_marker.ind])[0]
-    # subtract background from processed channel and if the background intensity for a certain pixel was larger than
-    # the processed channel, set intensity to 0 (no negative values)
-    back = np.rint(ne.evaluate("back * scalar"))
-    back = np.where(back>65535,65535,back.astype(np.uint16))
-    # set the pixel value to 0 if the image channel value is lower than the scaled background channel value
-    # otherwise, subtract.
-    output[channel] = np.where(image[channel]<back, 0, image[channel]-back)
-    back = None
-    return output[channel]
-
-def subtract(img, markers, output):
-    # iterating over channels
     for channel in range(len(markers)):
-        # this array is used to find the appropriate background channel
-        # IMPORTANT - all values in the 'marker_name' column should be unique        
-
-        # if no background channel is specified, channel is kept like if should be
         if markers.background.isnull()[channel] == True:
-            # no change - applies to background channels as well
-            output[channel] = copy.copy(img[channel])
             print(f"Channel {markers.marker_name[channel]} ({channel}) processed, no background subtraction")
-
         else:
-            # this array is used to find the appropriate background channel
-            # IMPORTANT - all values in the 'marker_name' column should be unique
-            find_background = np.array(markers.marker_name == markers.background[channel])
-            # scalar is the coefficient with which the background pixel intensity is scaled
-            # it is the result of dividing the marker exposure time with the background exposure time
-            # Marker_corr = Marker_root - Background * Exposure_marker / Exposure_background
-            background_marker = markers.iloc[find_background]
-
-            output[channel] = subtract_channel(img, markers, channel, background_marker, output)
-
-            print(f"Channel {markers.marker_name[channel]} ({channel}) processed, {markers.background[channel]} background channel subtracted")
-    return output
-
-def remove_back(img, markers):
-    #with dask.config.set(**{'array.slicing.split_large_chunks': True}):
-    # subset the image based on the remove column in the markers file
-    img = img[markers.remove.tolist()]
-    print()
-    print(f'Image shape: {img.shape}')
-    return img
-    
-def subres_tiles(level, level_full_shapes, tile_shapes, outpath, scale):
-    print(f"\n processing level {level}")
-    assert level >= 1
-    num_channels, h, w = level_full_shapes[level]
-    tshape = tile_shapes[level] or (h, w)
-    tiff = tifffile.TiffFile(outpath)
-    zimg = zarr.open(tiff.aszarr(series=0, level=level-1, squeeze=False))
-    for c in range(num_channels):
-        sys.stdout.write(
-            f"\r  processing channel {c + 1}/{num_channels}"
-        )
-        sys.stdout.flush()
-        th = tshape[0] * scale
-        tw = tshape[1] * scale
-        for y in range(0, zimg.shape[1], th):
-            for x in range(0, zimg.shape[2], tw):
-                a = zimg[c, y:y+th, x:x+tw, 0]
-                a = skimage.transform.downscale_local_mean(
-                    a, (scale, scale)
-                )
-                if np.issubdtype(zimg.dtype, np.integer):
-                    a = np.around(a)
-                a = a.astype('uint16')
-                yield a
-
-def main(args):
-    img_raw = AI.AICSImage(args.root)
-    img = img_raw.get_image_dask_data("CYX")
-
-    markers_raw = pd.read_csv(args.markers)
-    markers = process_markers(copy.copy(markers_raw))
-
-    output = dask.array.empty_like(img)
-
-    output = subtract(img, markers, output)
-    output = remove_back(output, markers)
-
-    markers_raw = markers_raw[markers_raw.remove != True]
-    markers_raw = markers_raw.drop("remove", axis = 1)
-    markers_raw.to_csv(args.markerout, index=False)
+            background_marker = markers.iloc[np.array(markers.marker_name == markers.background[channel])]
+            scalar = markers[markers.ind == channel].exposure.values / background_marker.exposure.values
+            mosaic_out[channel] = subtract_channels(mosaic[channel], mosaic[background_marker.ind.values[0]], scalar, chunk_size)
+            print(f"Channel {markers.marker_name[channel]} ({channel}) processed, background subtraction")
+
+    # removes channels from the image as specified in the markers file
+    mosaic_out = mosaic_out[np.where(markers.keep)[0]]
+    channel_names = list(markers.marker_name[markers.keep])
+
+    write_pyramid(
+        [mosaic_out], out_path, channel_names=channel_names, downscale_factor=2, pixel_size=pixel_size, save_RAM=False, tile_size=args.tile_size
+    )
+    markers = markers[markers.keep]
+    markers = markers.drop(columns=['keep','ind'])
+    markers.to_csv(args.markerout, index=False)
 
-    # Processing metadata - highly adapted to Lunaphore outputs
-    # check if metadata is present
-    try:
-        print(img_raw.metadata.images[0])
-        metadata = img_raw.metadata
-        metadata = process_metadata(img_raw.metadata, markers)
-    except:
-        metadata = None   
-
-    if args.pixel_size != None:
-        # If specified, the input pixel size is used
-        pixel_size = args.pixel_size   
-
-    else:
-        try:
-            if img_raw.metadata.images[0].pixels.physical_size_x != None:
-                pixel_size = img_raw.metadata.images[0].pixels.physical_size_x
-            else:
-                pixel_size = 1.0
-        except:
-            # If no pixel size specified anywhere, use default 1.0
-            pixel_size = 1.0
-    if (args.pyramid == True) and (int(args.tile_size)<=max(output[0].shape)):
-            
-        # construct levels
-        tile_size = int(args.tile_size)
-        scale = 2
-
-        print()
-        dtype = output.dtype
-        base_shape = output[0].shape
-        num_channels = output.shape[0]
-        num_levels = (np.ceil(np.log2(max(base_shape) / tile_size)) + 1).astype(int)
-        factors = 2 ** np.arange(num_levels)
-        shapes = (np.ceil(np.array(base_shape) / factors[:,None])).astype(int)
-        print(base_shape)
-        print(np.ceil(np.log2(max(base_shape)/tile_size))+1)
-
-        print("Pyramid level sizes: ")
-        for i, shape in enumerate(shapes):
-            print(f"   level {i+1}: {format_shape(shape)}", end="")
-            if i == 0:
-                print("(original size)", end="")
-            print()
-        print()
-        print(shapes)  
-
-        level_full_shapes = []
-        for shape in shapes:
-            level_full_shapes.append((num_channels, shape[0], shape[1]))
-        level_shapes = shapes
-        tip_level = np.argmax(np.all(level_shapes < tile_size, axis=1))
-        tile_shapes = [
-            (tile_size, tile_size) if i <= tip_level else None
-            for i in range(len(level_shapes))
-        ]
-        
-        # write pyramid
-        with tifffile.TiffWriter(args.output, ome=True, bigtiff=True) as tiff:
-            tiff.write(
-                data = output,
-                shape = level_full_shapes[0],
-                subifds=int(num_levels-1),
-                dtype=dtype,
-                resolution=(10000 / pixel_size, 10000 / pixel_size, "centimeter"),
-                tile=tile_shapes[0]
-            )
-            for level, (shape, tile_shape) in enumerate(
-                    zip(level_full_shapes[1:], tile_shapes[1:]), 1
-            ):
-                tiff.write(
-                    data = subres_tiles(level, level_full_shapes, tile_shapes, args.output, scale),
-                    shape=shape,
-                    subfiletype=1,
-                    dtype=dtype,
-                    tile=tile_shape
-                )
-    else:
-        # write image
-        with tifffile.TiffWriter(args.output, ome=True, bigtiff=True) as tiff:
-            tiff.write(
-                data = output,
-                shape = output.shape,
-                dtype=output.dtype,
-                resolution=(10000 / pixel_size, 10000 / pixel_size, "centimeter"),
-            )
-    try:
-        tifffile.tiffcomment(args.output, to_xml(metadata))
-    except:
-        pass
-    # note about metadata: the channels, planes etc were adjusted not to include the removed channels, however
-    # the channel ids have stayed the same as before removal. E.g if channels 1 and 2 are removed,
-    # the channel ids in the metadata will skip indices 1 and 2 (channel_id:0, channel_id:3, channel_id:4 ...)
-    print()
-
-        
-    
 if __name__ == '__main__':
-    # Read in arguments
+    _version = 'v0.4.1'
     args = get_args()
 
     # Run script
diff --git a/environment.yml b/environment.yml
index e17117e..957a760 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,14 +1,18 @@
-name: bsub
+name: backsub
 channels:
- - anaconda
- - conda-forge
- - defaults
+  - conda-forge
+  - defaults
+  - anaconda
 dependencies:
- - libgcc-ng
- - python=3.9.13
- - aicsimageio=4.9.2
- - scikit-image=0.19.2
- - numexpr=2.8.3
- - tifffile=2022.8.12
- - scipy=1.9.3
- - procps-ng
+  - "python=3.9"
+  - "openslide=3.4.1"
+  - "scikit-image=0.19.2" 
+  - "numexpr=2.8.3" 
+  - "tifffile=2022.8.12" 
+  - "scipy=1.9.3"
+  - "pandas=2.1.1" 
+  - "zarr=2.3.2"
+  - procps-ng
+  - pip
+  - pip:
+    - palom
\ No newline at end of file
diff --git a/scripts/old/background_sub_v034.py b/scripts/old/background_sub_v034.py
new file mode 100644
index 0000000..69bdd6d
--- /dev/null
+++ b/scripts/old/background_sub_v034.py
@@ -0,0 +1,313 @@
+from __future__ import print_function, division
+from multiprocessing.spawn import import_main_path
+import sys
+import copy
+import argparse
+import numpy as np
+import tifffile
+import zarr
+import skimage.transform
+from aicsimageio import aics_image as AI
+import pandas as pd
+import numexpr as ne
+from ome_types import from_tiff, to_xml
+from os.path import abspath
+from argparse import ArgumentParser as AP
+import time
+import dask
+from memory_profiler import profile
+# This API is apparently changing in skimage 1.0 but it's not clear to
+# me what the replacement will be, if any. We'll explicitly import
+# this so it will break loudly if someone tries this with skimage 1.0.
+try:
+    from skimage.util.dtype import _convert as dtype_convert
+except ImportError:
+    from skimage.util.dtype import convert as dtype_convert
+
+
+# arg parser
+def get_args():
+    # Script description
+    description="""Subtracts background - Lunaphore platform"""
+
+    # Add parser
+    parser = AP(description=description, formatter_class=argparse.RawDescriptionHelpFormatter)
+
+    # Sections
+    inputs = parser.add_argument_group(title="Required Input", description="Path to required input file")
+    inputs.add_argument("-r", "--root", dest="root", action="store", required=True, help="File path to root image file.")
+    inputs.add_argument("-m", "--markers", dest="markers", action="store", required=True, help="File path to required markers.csv file")
+    inputs.add_argument("--pixel-size", metavar="SIZE", dest = "pixel_size", type=float, default = None, action = "store",help="pixel size in microns; default is 1.0")
+    inputs.add_argument("--pyramid", dest="pyramid", required=False, default=True, help="Should output be pyramidal?")
+    inputs.add_argument("--tile-size", dest="tile_size", required=False, default=1024, help="Tile size for pyramid generation")
+    inputs.add_argument("--version", action="version", version="v0.3.4")
+    
+    outputs = parser.add_argument_group(title="Output", description="Path to output file")
+    outputs.add_argument("-o", "--output", dest="output", action="store", required=True, help="Path to output file")
+    outputs.add_argument("-mo", "--marker-output", dest="markerout", action="store", required=True, help="Path to output marker file")
+
+    arg = parser.parse_args()
+
+    # Standardize paths
+    arg.root = abspath(arg.root)
+    arg.markers = abspath(arg.markers)
+    arg.output = abspath(arg.output)
+    arg.markerout = abspath(arg.markerout)
+
+    return arg
+
+def preduce(coords, img_in, img_out):
+    print(img_in.dtype)
+    (iy1, ix1), (iy2, ix2) = coords
+    (oy1, ox1), (oy2, ox2) = np.array(coords) // 2
+    tile = skimage.img_as_float32(img_in[iy1:iy2, ix1:ix2])
+    tile = skimage.transform.downscale_local_mean(tile, (2, 2))
+    tile = dtype_convert(tile, 'uint16')
+    #tile = dtype_convert(tile, img_in.dtype)
+    img_out[oy1:oy2, ox1:ox2] = tile
+
+def format_shape(shape):
+    return "%dx%d" % (shape[1], shape[0])
+
+def process_markers(markers):
+    # add column with starting indices (which match the image channel indices)
+    # this should be removed soon
+    markers['ind'] = range(0, len(markers))
+
+    # if the 'remove' column is not specified, all channels are kept. If it is 
+    # present, it is converted to a boolean indicating which channels should be removed
+    if 'remove' not in markers:
+        markers['remove'] = ["False" for i in range(len(markers))]
+    else:
+        markers['remove'] = markers['remove'] == True
+    # invert the markers.remove column to indicate which columns to keep
+    markers['remove'] = markers['remove'] == False
+
+    return markers
+
+def process_metadata(metadata, markers):
+    try:
+        metadata.screens[0].reagents = [metadata.screens[0].reagents[i] for i in markers[markers.remove == True].ind]
+    except:
+        pass
+    try:
+        metadata.structured_annotations = [metadata.structured_annotations[i] for i in markers[markers.remove == True].ind]
+    except:
+        pass
+    # these two are required
+    metadata.images[0].pixels.size_c = len(markers[markers.remove == True])
+    metadata.images[0].pixels.channels = [metadata.images[0].pixels.channels[i] for i in markers[markers.remove == True].ind]
+    try:
+        metadata.images[0].pixels.planes = [metadata.images[0].pixels.planes[i] for i in markers[markers.remove == True].ind]
+    except:
+        pass
+    try:
+        metadata.images[0].pixels.tiff_data_blocks[0].plane_count = sum(markers.remove == True)
+    except:
+        pass
+    return metadata
+
+
+# NaN values return True for the statement below in this version of Python. Did not use math.isnan() since the values
+# are strings if present
+def isNaN(x):
+    return x != x
+
+def subtract_channel(image, markers, channel, background_marker, output):
+    scalar = markers[markers.ind == channel].exposure.values / background_marker.exposure.values
+    
+    # create temporary dataframe which will store the multiplied background rounded up to nearest integer
+    # [0] at the end needed to get [x, y] shape, and not have [1, x, y]
+    back = copy.copy(image[background_marker.ind])[0]
+    # subtract background from processed channel and if the background intensity for a certain pixel was larger than
+    # the processed channel, set intensity to 0 (no negative values)
+    back = np.rint(ne.evaluate("back * scalar"))
+    back = np.where(back>65535,65535,back.astype(np.uint16))
+    # set the pixel value to 0 if the image channel value is lower than the scaled background channel value
+    # otherwise, subtract.
+    output[channel] = np.where(image[channel]<back, 0, image[channel]-back)
+    back = None
+    return output[channel]
+
+def subtract(img, markers, output):
+    # iterating over channels
+    for channel in range(len(markers)):
+        # this array is used to find the appropriate background channel
+        # IMPORTANT - all values in the 'marker_name' column should be unique        
+
+        # if no background channel is specified, channel is kept like if should be
+        if markers.background.isnull()[channel] == True:
+            # no change - applies to background channels as well
+            output[channel] = copy.copy(img[channel])
+            print(f"Channel {markers.marker_name[channel]} ({channel}) processed, no background subtraction")
+
+        else:
+            # this array is used to find the appropriate background channel
+            # IMPORTANT - all values in the 'marker_name' column should be unique
+            find_background = np.array(markers.marker_name == markers.background[channel])
+            # scalar is the coefficient with which the background pixel intensity is scaled
+            # it is the result of dividing the marker exposure time with the background exposure time
+            # Marker_corr = Marker_root - Background * Exposure_marker / Exposure_background
+            background_marker = markers.iloc[find_background]
+
+            output[channel] = subtract_channel(img, markers, channel, background_marker, output)
+
+            print(f"Channel {markers.marker_name[channel]} ({channel}) processed, {markers.background[channel]} background channel subtracted")
+    return output
+
+def remove_back(img, markers):
+    #with dask.config.set(**{'array.slicing.split_large_chunks': True}):
+    # subset the image based on the remove column in the markers file
+    img = img[markers.remove.tolist()]
+    print()
+    print(f'Image shape: {img.shape}')
+    return img
+    
+def subres_tiles(level, level_full_shapes, tile_shapes, outpath, scale):
+    print(f"\n processing level {level}")
+    assert level >= 1
+    num_channels, h, w = level_full_shapes[level]
+    tshape = tile_shapes[level] or (h, w)
+    tiff = tifffile.TiffFile(outpath)
+    zimg = zarr.open(tiff.aszarr(series=0, level=level-1, squeeze=False))
+    for c in range(num_channels):
+        sys.stdout.write(
+            f"\r  processing channel {c + 1}/{num_channels}"
+        )
+        sys.stdout.flush()
+        th = tshape[0] * scale
+        tw = tshape[1] * scale
+        for y in range(0, zimg.shape[1], th):
+            for x in range(0, zimg.shape[2], tw):
+                a = zimg[c, y:y+th, x:x+tw, 0]
+                a = skimage.transform.downscale_local_mean(
+                    a, (scale, scale)
+                )
+                if np.issubdtype(zimg.dtype, np.integer):
+                    a = np.around(a)
+                a = a.astype('uint16')
+                yield a
+
+@profile
+def main(args):
+    img_raw = AI.AICSImage(args.root)
+    img = img_raw.get_image_dask_data("CYX")
+
+    markers_raw = pd.read_csv(args.markers)
+    markers = process_markers(copy.copy(markers_raw))
+
+    output = dask.array.empty_like(img)
+
+    output = subtract(img, markers, output)
+    output = remove_back(output, markers)
+
+    markers_raw = markers_raw[markers_raw.remove != True]
+    markers_raw = markers_raw.drop("remove", axis = 1)
+    markers_raw.to_csv(args.markerout, index=False)
+
+    # Processing metadata - highly adapted to Lunaphore outputs
+    # check if metadata is present
+    try:
+        print(img_raw.metadata.images[0])
+        metadata = img_raw.metadata
+        metadata = process_metadata(img_raw.metadata, markers)
+    except:
+        metadata = None   
+
+    if args.pixel_size != None:
+        # If specified, the input pixel size is used
+        pixel_size = args.pixel_size   
+
+    else:
+        try:
+            if img_raw.metadata.images[0].pixels.physical_size_x != None:
+                pixel_size = img_raw.metadata.images[0].pixels.physical_size_x
+            else:
+                pixel_size = 1.0
+        except:
+            # If no pixel size specified anywhere, use default 1.0
+            pixel_size = 1.0
+    if (args.pyramid == True) and (int(args.tile_size)<=max(output[0].shape)):
+            
+        # construct levels
+        tile_size = int(args.tile_size)
+        scale = 2
+
+        print()
+        dtype = output.dtype
+        base_shape = output[0].shape
+        num_channels = output.shape[0]
+        num_levels = (np.ceil(np.log2(max(base_shape) / tile_size)) + 1).astype(int)
+        factors = 2 ** np.arange(num_levels)
+        shapes = (np.ceil(np.array(base_shape) / factors[:,None])).astype(int)
+        print(base_shape)
+        print(np.ceil(np.log2(max(base_shape)/tile_size))+1)
+
+        print("Pyramid level sizes: ")
+        for i, shape in enumerate(shapes):
+            print(f"   level {i+1}: {format_shape(shape)}", end="")
+            if i == 0:
+                print("(original size)", end="")
+            print()
+        print()
+        print(shapes)  
+
+        level_full_shapes = []
+        for shape in shapes:
+            level_full_shapes.append((num_channels, shape[0], shape[1]))
+        level_shapes = shapes
+        tip_level = np.argmax(np.all(level_shapes < tile_size, axis=1))
+        tile_shapes = [
+            (tile_size, tile_size) if i <= tip_level else None
+            for i in range(len(level_shapes))
+        ]
+        
+        # write pyramid
+        with tifffile.TiffWriter(args.output, ome=True, bigtiff=True) as tiff:
+            tiff.write(
+                data = output,
+                shape = level_full_shapes[0],
+                subifds=int(num_levels-1),
+                dtype=dtype,
+                resolution=(10000 / pixel_size, 10000 / pixel_size, "centimeter"),
+                tile=tile_shapes[0]
+            )
+            for level, (shape, tile_shape) in enumerate(
+                    zip(level_full_shapes[1:], tile_shapes[1:]), 1
+            ):
+                tiff.write(
+                    data = subres_tiles(level, level_full_shapes, tile_shapes, args.output, scale),
+                    shape=shape,
+                    subfiletype=1,
+                    dtype=dtype,
+                    tile=tile_shape
+                )
+    else:
+        # write image
+        with tifffile.TiffWriter(args.output, ome=True, bigtiff=True) as tiff:
+            tiff.write(
+                data = output,
+                shape = output.shape,
+                dtype=output.dtype,
+                resolution=(10000 / pixel_size, 10000 / pixel_size, "centimeter"),
+            )
+    try:
+        tifffile.tiffcomment(args.output, to_xml(metadata))
+    except:
+        pass
+    # note about metadata: the channels, planes etc were adjusted not to include the removed channels, however
+    # the channel ids have stayed the same as before removal. E.g if channels 1 and 2 are removed,
+    # the channel ids in the metadata will skip indices 1 and 2 (channel_id:0, channel_id:3, channel_id:4 ...)
+    print()
+
+        
+    
+if __name__ == '__main__':
+    # Read in arguments
+    args = get_args()
+
+    # Run script
+    st = time.time()
+    main(args)
+    rt = time.time() - st
+    print(f"Script finished in {rt // 60:.0f}m {rt % 60:.0f}s")