From 330abf8bbc07c29be736033bb4baed79fe1dc438 Mon Sep 17 00:00:00 2001
From: Christopher Rhodes <christopher.rhodes@embl.de>
Date: Tue, 26 Sep 2023 10:29:40 +0200
Subject: [PATCH] Roughed in chaetoceros project
---
extensions/{examples => chaeo}/__init__.py | 0
extensions/chaeo/products.py | 97 ++++++++++++++++
extensions/chaeo/test_zstack.py | 31 +++++
extensions/chaeo/zstack.py | 109 ++++++++++++++++++
extensions/ilastik/examples/__init__.py | 0
.../{ => ilastik}/examples/ilastik3d.py | 0
extensions/ilastik/validators.py | 5 +
model_server/batch.py | 22 ++++
model_server/process.py | 52 +++++++++
readme.md | 17 +++
10 files changed, 333 insertions(+)
rename extensions/{examples => chaeo}/__init__.py (100%)
create mode 100644 extensions/chaeo/products.py
create mode 100644 extensions/chaeo/test_zstack.py
create mode 100644 extensions/chaeo/zstack.py
create mode 100644 extensions/ilastik/examples/__init__.py
rename extensions/{ => ilastik}/examples/ilastik3d.py (100%)
create mode 100644 extensions/ilastik/validators.py
create mode 100644 model_server/batch.py
create mode 100644 model_server/process.py
create mode 100644 readme.md
diff --git a/extensions/examples/__init__.py b/extensions/chaeo/__init__.py
similarity index 100%
rename from extensions/examples/__init__.py
rename to extensions/chaeo/__init__.py
diff --git a/extensions/chaeo/products.py b/extensions/chaeo/products.py
new file mode 100644
index 00000000..a0f50de3
--- /dev/null
+++ b/extensions/chaeo/products.py
@@ -0,0 +1,97 @@
+from PIL import Image, ImageDraw, ImageFont
+
+def draw_boxes_on_2d_image(img, boxes, **kwargs):
+ pilimg = Image.fromarray(np.copy(img)) # drawing modifies array in-place
+ draw = ImageDraw.Draw(pilimg)
+ font_size = kwargs.get('font_size', 18)
+ linewidth = kwargs.get('linewidth', 4)
+
+ draw.font = ImageFont.truetype(font="arial.ttf", size=font_size)
+
+ for box in boxes:
+ y0 = box['info'].y0
+ y1 = box['info'].y1
+ x0 = box['info'].x0
+ x1 = box['info'].x1
+ xm = round((x0 + x1) / 2)
+
+ la = box['info'].label
+ zi = box['info'].zi
+
+ draw.rectangle([(x0, y0), (x1, y1)], outline='white', width=linewidth)
+
+ if kwargs.get('add_label') is True:
+ draw.text((xm, y0), f'Z{zi:04d}-L{la:04d}', fill='white', anchor='mb')
+
+ return pilimg
+
+def generate_patches(
+ desc, stack, boxes, rescale_clip=0.0,
+ pad_to=256,
+ proj=lambda x: x.max(axis=0),
+ prefix='patch',
+ **kwargs
+):
+ patch_dir = root / 'output' / 'patches' / desc
+ patch_dir.mkdir(parents=True, exist_ok=True)
+
+ for box in boxes:
+ obj = box['info']
+ sl = box['slice']
+ rbb = box['relative_bounding_box']
+
+ patch = proj(stack[sl])
+ patch_fname = f'{prefix}-la{obj.label:04d}-zi{obj.zi:04d}'
+
+ if rescale_clip is not None:
+ patch = rescale(patch, rescale_clip)
+
+ if kwargs.get('draw_bounding_box') is True:
+ x0 = rbb['x0']
+ y0 = rbb['y0']
+ x1 = rbb['x1']
+ y1 = rbb['y1']
+
+ pilimg = Image.fromarray(patch) # drawing modifies array in-place
+ draw = ImageDraw.Draw(pilimg)
+ draw.rectangle([(x0, y0), (x1, y1)], outline='white', width=kwargs.get('linewidth', 1))
+ patch = np.array(pilimg)
+
+ if pad_to:
+ patch = pad(patch, pad_to)
+
+ imsave(
+ patch_dir / (patch_fname + '.png'),
+ resample(patch),
+ check_contrast=False,
+ )
+ print(f'Successfully wrote {len(boxes)} patches to:\n{patch_dir}')
+
+
+def generate_3d_patches( # in extensions.chaeo.products
+ desc, stack, boxes, rescale_clip=0.0,
+ pad_to=256,
+ prefix='patch',
+ proj=lambda x: x,
+):
+ patch_dir = root / 'output' / '3d_patches' / desc
+ patch_dir.mkdir(parents=True, exist_ok=True)
+
+ for box in boxes:
+ obj = box['info']
+ sl = box['slice']
+ patch = proj(stack[sl])
+ patch_fname = f'{prefix}-la{obj.label:04d}-zi{obj.zi:04d}'
+
+ if rescale_clip is not None:
+ patch = rescale(patch, rescale_clip)
+
+ if pad_to:
+ patch = pad_3d(patch, pad_to)
+
+ imwrite(
+ patch_dir / (patch_fname + '.tif'),
+ patch,
+ imagej=True
+ )
+ print(f'Successfully wrote {len(boxes)} patches to:\n{patch_dir}')
\ No newline at end of file
diff --git a/extensions/chaeo/test_zstack.py b/extensions/chaeo/test_zstack.py
new file mode 100644
index 00000000..1f9fcf96
--- /dev/null
+++ b/extensions/chaeo/test_zstack.py
@@ -0,0 +1,31 @@
+import unittest
+
+from extensions.chaeo.zstack import build_stack_mask
+from model_server.accessors import InMemoryDataAccessor
+
+class TestZStackDerivedDataProducts(unittest.TestCase):
+
+ def setUp(self) -> None:
+ # need test data incl obj map
+ self.obmap = None
+ self.stack = None
+
+ def test_zmask_makes_correct_boxes(self):
+ zmask, meta = build_stack_mask(
+ 'test_zmask_with boxes',
+ self.stack,
+ self.obmap,
+ mask_type='boxes',
+ )
+ zmask_acc = InMemoryDataAccessor(zmask)
+ self.assertTrue(zmask_acc.is_object_map())
+
+ # assert dimensionality of zmask
+ self.assertEqual(zmask.shape_dict['Z'] > 1)
+ self.assertEqual(zmask.shape_dict['C'] == 1)
+
+ # assert non-trivial meta info in boxes
+ pass
+
+ def test_zmask_makes_correct_contours(self):
+ pass
\ No newline at end of file
diff --git a/extensions/chaeo/zstack.py b/extensions/chaeo/zstack.py
new file mode 100644
index 00000000..475f6fd8
--- /dev/null
+++ b/extensions/chaeo/zstack.py
@@ -0,0 +1,109 @@
+import numpy as np
+import pandas as pd
+
+from skimage.measure import find_contours, regionprops_table
+
+# build a single boolean 3d mask (objects v. bboxes) and return bounding boxes
+def build_stack_mask(desc, obmap, stack, filters=None, mask_type='contour', expand_box_by=(0, 0)): # TODO: specify boxes data type
+ """
+
+ filters: dict of (min, max) tuples
+ expand_box_by: (xy, z) pixelsf
+ """
+
+ # validate inputs
+ assert len(stack.shape) == 3, stack.shape
+ assert mask_type in ('contour', 'box'), mask_type # TODO: replace with call to validator
+
+ for k in filters.keys():
+ assert k in ('area', 'solidity')
+ vmin, vmax = filters[k]
+ assert vmin >= 0
+
+ # build object query
+ query_str = 'label > 0' # always true
+ if filters:
+ for k in filters.keys():
+ assert k in ('area', 'solidity')
+ vmin, vmax = filters[k]
+ assert vmin >= 0
+ query_str = query_str + f' & {k} > {vmin} & {k} < {vmax}'
+
+ # build dataframe of objects, assign z index to each object
+ argmax = stack.argmax(axis=0)
+ df = (
+ pd.DataFrame(
+ regionprops_table(
+ obmap,
+ intensity_image=argmax,
+ properties=('label', 'area', 'intensity_mean', 'solidity', 'bbox')
+ )
+ )
+ .query(query_str)
+ .rename(
+ columns={
+ 'bbox-0': 'y0',
+ 'bbox-1': 'x0',
+ 'bbox-2': 'y1',
+ 'bbox-3': 'x1',
+ }
+ )
+ )
+ df['zi'] = df['intensity_mean'].round().astype('int')
+
+ # make an object map where label is replaced by focus position in stack and background is -1
+ lut = np.zeros(obmap.max() + 1) - 1
+ lut[df.label] = df.zi
+
+ # convert bounding boxes to slices
+ ebxy, ebz = expand_box_by
+ nz, h, w = stack.shape
+
+ boxes = []
+ for ob in df.itertuples(name='LabeledObject'):
+ y0 = max(ob.y0 - ebxy, 0)
+ y1 = min(ob.y1 + ebxy, h - 1)
+ x0 = max(ob.x0 - ebxy, 0)
+ x1 = min(ob.x1 + ebxy, w - 1)
+ z0 = max(ob.zi - ebz, 0)
+ z1 = min(ob.zi + ebz, nz)
+
+ # relative bounding box positions
+ rbb = {
+ 'y0': ob.y0 - y0,
+ 'y1': ob.y1 - y0,
+ 'x0': ob.x0 - x0,
+ 'x1': ob.x1 - x0,
+ }
+
+ sl = np.s_[z0: z1 + 1, y0: y1, x0: x1]
+
+ # compute contours
+ obmask = (obmap == ob.label)
+ contour = find_contours(obmask)
+ mask = obmask[ob.y0: ob.y1, ob.x0: ob.x1]
+
+ boxes.append({
+ 'info': ob,
+ 'slice': sl,
+ 'relative_bounding_box': rbb,
+ 'contour': contour,
+ 'mask': mask
+ })
+
+ # build mask z-stack
+ zi_st = np.zeros(stack.shape, dtype='bool')
+ if mask_type == 'contour':
+ zi_map = (lut[obmap] + 1.0).astype('int')
+ idxs = np.array([zi_map]) - 1
+ np.put_along_axis(zi_st, idxs, 1, axis=0)
+
+ # change background level from to 0 in final frame
+ zi_st[-1, :, :][obmap == 0] = 0
+
+ elif mask_type == 'box':
+ for bb in boxes:
+ sl = bb['slice']
+ zi_st[sl] = 1
+
+ return zi_st, boxes
\ No newline at end of file
diff --git a/extensions/ilastik/examples/__init__.py b/extensions/ilastik/examples/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/extensions/examples/ilastik3d.py b/extensions/ilastik/examples/ilastik3d.py
similarity index 100%
rename from extensions/examples/ilastik3d.py
rename to extensions/ilastik/examples/ilastik3d.py
diff --git a/extensions/ilastik/validators.py b/extensions/ilastik/validators.py
new file mode 100644
index 00000000..d82e24a8
--- /dev/null
+++ b/extensions/ilastik/validators.py
@@ -0,0 +1,5 @@
+from model_server.accessors import GenericImageDataAccessor
+
+def is_object_map(img: GenericImageDataAccessor):
+ # TODO: implement
+ pass
\ No newline at end of file
diff --git a/model_server/batch.py b/model_server/batch.py
new file mode 100644
index 00000000..999aa4cc
--- /dev/null
+++ b/model_server/batch.py
@@ -0,0 +1,22 @@
+def get_tif_seq(where, prefix=None, verbose=True):
+ files = [ff for ff in os.listdir(where) if ff.startswith(prefix)]
+ files.sort()
+ paths = [where / fn for fn in files]
+ ndas = [imread(pa.__str__()) for pa in paths]
+ assert all([n.shape == ndas[0].shape and n.dtype == ndas[0].dtype for n in ndas])
+ if verbose:
+ print(f'Successfully read {len(ndas)} x {ndas[0].shape} px images of type {ndas[0].dtype} from:\n{where}')
+ return ndas
+
+def write_tif_zstack(data, desc, dtype='uint16', subdir='output'):
+ where = root / subdir
+ filename = desc + '.tif'
+ if not os.path.exists(where):
+ os.makedirs(where)
+ abspath = where / filename
+ imwrite(
+ abspath,
+ np.array(data).astype(dtype),
+ imagej=True
+ )
+ print(f'Successfully wrote {len(data)} x {data[0].shape} px images of type {dtype} to:\n{abspath}')
\ No newline at end of file
diff --git a/model_server/process.py b/model_server/process.py
new file mode 100644
index 00000000..26fd3005
--- /dev/null
+++ b/model_server/process.py
@@ -0,0 +1,52 @@
+"""
+Image processing utility functions
+"""
+from math import ceil, floor
+
+import numpy as np
+from skimage.exposure import rescale_intensity
+
+
+def pad(im, mpx): # now in model_server.batch
+ '''Pads and crops image width edge values to specified dimension'''
+ dh = 0.5 * (mpx - im.shape[0])
+ dw = 0.5 * (mpx - im.shape[1])
+
+ if dw < 0:
+ x0 = floor(-dw)
+ x1 = x0 + mpx
+ im = im[:, x0:x1]
+ dw = 0
+ if dh < 0:
+ y0 = floor(-dh)
+ y1 = y0 + mpx
+ im = im[y0:y1, :]
+ dh = 0
+
+ border = ((floor(dh), ceil(dh)), (floor(dw), ceil(dw)))
+ padded = np.pad(im, border, mode='constant')
+ if padded.shape != (mpx, mpx):
+ raise Exception(f'Incorrect image shape: {padded.shape} v. {(mpx, mpx)}')
+ return padded
+
+def pad_3d(im, mpx): # im: [z x h x w]
+ assert(len(im.shape) == 3)
+ nz, h, w = im.shape
+ padded = np.zeros((nz, mpx, mpx), dtype=im.dtype)
+ for zi in range(nz):
+ padded[zi, :, :] = pad(im[zi, :, :], mpx)
+ return padded
+
+def resample(nda, cmin=0, cmax=2**16): # now in model_server.batch
+ return rescale_intensity(
+ np.clip(nda, cmin, cmax),
+ in_range=(cmin, cmax + 1),
+ out_range=(0, 2**8)
+ ).astype('uint8')
+
+
+def rescale(nda, clip=0.0): # now in model_server.batch
+ clip_pct = (100.0 * clip, 100.0 * (1.0 - clip))
+ cmin, cmax = np.percentile(nda, clip_pct)
+ rescaled = rescale_intensity(nda, in_range=(cmin, cmax))
+ return rescaled
\ No newline at end of file
diff --git a/readme.md b/readme.md
new file mode 100644
index 00000000..c424e968
--- /dev/null
+++ b/readme.md
@@ -0,0 +1,17 @@
+# model_server
+
+# How to extend service
+Add sub-package to extensions
+Add models that inherit from model_server.Model
+In workflows, implement pipelines with File I/O via accessors.GenericImageDataAccessor
+(to decouple model logic from image data source)
+Set extensions-specific folders, etc. in conf relative to overall package root (set by user)
+As much as possible, set pipeline and model parameters with defaults and support overrides by optional API arguments;
+this helps non-coding users control their jobs
+Set up API endpoints in router, following as much as possible existing conventions with load, infer, etc. keyword
+
+decouple data access from processing
+
+control either via batch runners or API (serial)
+
+workflow: combines data access with processing via models, produces primary outputs
\ No newline at end of file
--
GitLab