Batch-export 3d patches and table of focus metrics; optionally annotate...

Batch-export 3d patches and table of focus metrics; optionally annotate bounding box at plane of best focus

extensions/chaeo/examples/batch_run_patches.py

-import json
 from pathlib import Path
 import re
 from time import localtime, strftime
@@ -10,7 +9,7 @@ from model_server.accessors import InMemoryDataAccessor, write_accessor_data_to_
 
 if __name__ == '__main__':
     where_czi = Path(
-        'z:/rhodes/projects/proj0004-marine-photoactivation/data/exp0038/AutoMic/20230906-163415/Selection'
+        'c:/Users/rhodes/projects/proj0004-marine-photoactivation/data/exp0038/AutoMic/20230906-163415/Selection'
     )
 
     where_output_root = Path(
@@ -67,5 +66,4 @@ if __name__ == '__main__':
             write_accessor_data_to_file(
                 where_output / k / (ff.stem + '.tif'),
                 InMemoryDataAccessor(result['interm'][k])
-            )
-        break
\ No newline at end of file
+            )
\ No newline at end of file

extensions/chaeo/examples/export_patch_focus_metrics.py

+from pathlib import Path
+import re
+from time import localtime, strftime
+from typing import Dict
+
+import pandas as pd
+
+from extensions.ilastik.models import IlastikPixelClassifierModel
+from extensions.chaeo.products import export_3d_patches_with_focus_metrics
+from extensions.chaeo.zmask import build_zmask_from_object_mask
+from model_server.accessors import generate_file_accessor, InMemoryDataAccessor, write_accessor_data_to_file
+from model_server.workflows import Timer
+
+
+def export_patch_focus_metrics_from_multichannel_zstack(
+        input_zstack_path: str,
+        ilastik_project_file: str,
+        pxmap_threshold: float,
+        pixel_class: int,
+        zmask_channel: int,
+        patches_channel: int,
+        where_output: str,
+        mask_type: str = 'boxes',
+        zmask_filters: Dict = None,
+        zmask_expand_box_by: int = None,
+        annotate_focus_metric=None,
+        **kwargs,
+) -> Dict:
+
+    ti = Timer()
+    stack = generate_file_accessor(Path(input_zstack_path))
+    fstem = Path(input_zstack_path).stem
+    ti.click('file_input')
+    assert stack.nz > 1, 'Expecting z-stack'
+
+    # MIP and classify pixels
+    mip = InMemoryDataAccessor(
+        stack.get_one_channel_data(channel=0).data.max(axis=-1, keepdims=True)
+    )
+    px_model = IlastikPixelClassifierModel(
+        params={'project_file': Path(ilastik_project_file)}
+    )
+    pxmap, _ = px_model.infer(mip)
+    ti.click('infer_pixel_probability')
+
+    obmask = InMemoryDataAccessor(
+        pxmap.data > pxmap_threshold
+    )
+    ti.click('threshold_pixel_mask')
+
+    # make zmask
+    zmask, zmask_meta, df, interm = build_zmask_from_object_mask(
+        obmask.get_one_channel_data(pixel_class),
+        stack.get_one_channel_data(zmask_channel),
+        mask_type=mask_type,
+        filters=zmask_filters,
+        expand_box_by=zmask_expand_box_by,
+    )
+    zmask_acc = InMemoryDataAccessor(zmask)
+    ti.click('generate_zmasks')
+
+    files = export_3d_patches_with_focus_metrics(
+        Path(where_output) / '3d_patches',
+        stack.get_one_channel_data(patches_channel),
+        zmask_meta,
+        prefix=fstem,
+        draw_bounding_box=False,
+        rescale_clip=0.0,
+        make_3d=True,
+        annotate_focus_metric=annotate_focus_metric,
+    )
+    ti.click('export_patches')
+
+    return {
+        'pixel_model_id': px_model.model_id,
+        'input_filepath': input_zstack_path,
+        'number_of_objects': len(zmask_meta),
+        'success': True,
+        'timer_results': ti.events,
+        'dataframe': df,
+        'interm': interm,
+    }
+
+if __name__ == '__main__':
+    where_czi = Path(
+        'c:/Users/rhodes/projects/proj0004-marine-photoactivation/data/exp0038/AutoMic/20230906-163415/Selection'
+    )
+
+    where_output_root = Path(
+        'c:/Users/rhodes/projects/proj0011-plankton-seg/exp0009/output'
+    )
+    yyyymmdd = strftime('%Y%m%d', localtime())
+    idx = 0
+    while Path(where_output_root / f'batch-output-{yyyymmdd}-{idx:04d}').exists():
+        idx += 1
+    where_output = Path(
+        where_output_root / f'batch-output-{yyyymmdd}-{idx:04d}'
+    )
+
+    csv_args = {'mode': 'w', 'header': True} # when creating file
+    px_ilp = Path.home() / 'model-server' / 'ilastik' / 'AF405-bodies_boundaries.ilp'
+
+    for ff in where_czi.iterdir():
+        pattern = 'Selection--W([\d]+)--P([\d]+)-T([\d]+)'
+        ma = re.match(pattern, ff.stem)
+
+        print(ff)
+        if not ff.suffix.upper() == '.CZI':
+            continue
+        if int(ma.groups()[1]) > 10: # skip second half of set
+            continue
+
+        export_kwargs = {
+            'input_zstack_path': (where_czi / ff).__str__(),
+            'ilastik_project_file': px_ilp.__str__(),
+            'pxmap_threshold': 0.25,
+            'pixel_class': 0,
+            'zmask_channel': 0,
+            'patches_channel': 4,
+            'where_output': where_output.__str__(),
+            'mask_type': 'boxes',
+            'zmask_filters': {'area': (1e3, 1e8)},
+            'zmask_expand_box_by': (128, 3),
+            'annotate_focus_metric': 'rms_sobel'
+        }
+
+        result = export_patch_focus_metrics_from_multichannel_zstack(**export_kwargs)
+
+        # parse and record results
+        df = result['dataframe']
+        df['filename'] = ff.name
+        df.to_csv(where_output / 'df_objects.csv', **csv_args)
+        pd.DataFrame(result['timer_results'], index=[0]).to_csv(where_output / 'timer_results.csv', **csv_args)
+        pd.json_normalize(export_kwargs).to_csv(where_output / 'workflow_params.csv', **csv_args)
+        csv_args = {'mode': 'a', 'header': False} # append to CSV from here on
+
+        # export intermediate data if flagged
+        for k in result['interm'].keys():
+            write_accessor_data_to_file(
+                where_output / k / (ff.stem + '.tif'),
+                InMemoryDataAccessor(result['interm'][k])
+            )
\ No newline at end of file

extensions/chaeo/products.py

+from math import sqrt
 from pathlib import Path
 
 import numpy as np
-
+import pandas as pd
+from scipy.stats import moment
+from skimage.filters import sobel
 from skimage.io import imsave
+from skimage.measure import shannon_entropy
 from tifffile import imwrite
 
 from extensions.chaeo.annotators import draw_box_on_patch
@@ -80,3 +84,86 @@ def export_patches_from_zstack(
         _write_patch_to_file(where, fname, resample_to_8bit(patch))
         exported.append(fname)
     return exported
+
+def export_3d_patches_with_focus_metrics(
+        where: Path,
+        stack: GenericImageDataAccessor,
+        zmask_meta: list,
+        rescale_clip: float = 0.0,
+        pad_to: int = 256,
+        prefix='patch',
+        **kwargs
+):
+    """
+    Export 3D patches as multi-level z-stacks, along with CSV of various focus methods for each z-position
+
+    :param kwargs:
+        annotate_focus_metric: name focus metric to use when drawing bounding box at optimal focus z-position
+    :return:
+        list of exported files
+    """
+    assert stack.chroma == 1, 'Expecting monochromatic image data'
+    assert stack.nz > 1, 'Expecting z-stack'
+
+    def get_zstack_focus_metrics(zs):
+        nz = zs.shape[3]
+        me = {
+            'max_intensity': lambda x: np.max(x),
+            'stdev': lambda x: np.std(x),
+            'max_sobel': lambda x: np.max(sobel(x)),
+            'rms_sobel': lambda x: sqrt(np.mean(sobel(x) ** 2)),
+            'entropy': lambda x: shannon_entropy(x),
+            'moment': lambda x: moment(x, moment=2),
+        }
+        dd = {}
+        for zi in range(0, nz):
+            spf = zs[:, :, :, zi].flatten()
+            dd[zi] = {k: me[k](spf) for k in me.keys()}
+        return dd
+
+    exported = []
+    patch_meta = []
+    for mi in zmask_meta:
+        obj = mi['info']
+        sl = mi['slice']
+        rbb = mi['relative_bounding_box']
+
+        patch = stack.data[sl]
+
+        assert len(patch.shape) == 4
+        assert patch.shape[2] == stack.chroma
+
+        if rescale_clip is not None:
+            patch = rescale(patch, rescale_clip)
+
+        # unpack relative bounding box and define subset of patch data
+        x0 = rbb['x0']
+        y0 = rbb['y0']
+        x1 = rbb['x1']
+        y1 = rbb['y1']
+        sp_sl = np.s_[y0: y1, x0: x1, :, :]
+        subpatch = patch[sp_sl]
+
+        # compute focus metrics for all z-levels
+        me_dict = get_zstack_focus_metrics(subpatch)
+        patch_meta.append({'label': obj.label, 'zi': obj.zi, 'metrics': me_dict})
+        me_df = pd.DataFrame(me_dict).T
+
+        # drawing bounding box only on focused slice
+        if ak := kwargs.get('annotate_focus_metric') in me_dict.keys():
+            zi_foc = me_df.idxmax().to_dict()[ak]
+            patch[:, :, 0, zi_foc] = draw_box_on_patch(
+                patch[:, :, 0, zi_foc],
+                ((x0, y0), (x1, y1)),
+            )
+
+        if pad_to:
+            patch = pad(patch, pad_to)
+
+        fstem = f'{prefix}-la{obj.label:04d}-zi{obj.zi:04d}'
+        _write_patch_to_file(where, fstem + '.tif', resample_to_8bit(patch))
+        exported.append(fstem + '.tif')
+        me_df.to_csv(where / (fstem + '.csv'))
+        exported.append(fstem + '.csv')
+
+    return exported

extensions/chaeo/workflows.py

@@ -64,20 +64,26 @@ def export_patches_from_multichannel_zstack(
     import numpy as np
     from skimage.filters import gaussian, sobel
 
-    def zs_projector(zs):
-        sigma = 1.5
-        blur = gaussian(sobel(zs), sigma)
-        argmax = np.argmax(blur, axis=3, keepdims=True)
-        return np.take_along_axis(zs, argmax, axis=3)
+    # def zs_projector(zs):
+    #     sigma = 1.5
+    #     blur = gaussian(sobel(zs), sigma)
+    #     argmax = np.argmax(blur, axis=3, keepdims=True)
+    #     return np.take_along_axis(zs, argmax, axis=3)
+    #
+    # def zs_annotate_best_focus(zs):
+    #     pass
 
     files = export_patches_from_zstack(
-        Path(where_output) / '2d_patches',
+        # Path(where_output) / '2d_patches',
+        Path(where_output) / '3d_patches',
         stack.get_one_channel_data(patches_channel),
         zmask_meta,
         prefix=fstem,
-        draw_bounding_box=True,
+        # draw_bounding_box=True,
+        draw_bounding_box=False,
         rescale_clip=0.0,
-        projector=zs_projector,
+        # projector=zs_projector,
+        make_3d=True,
     )
     ti.click('export_patches')
 
@@ -94,15 +100,15 @@ def export_patches_from_multichannel_zstack(
     )
     ti.click('export_annotated_zstack')
 
-    # generate multichannel projection from label centroids
-    dff = df[df['keeper']]
-    interm['projected'] = project_stack_from_focal_points(
-        dff['centroid-0'].to_numpy(),
-        dff['centroid-1'].to_numpy(),
-        dff['zi'].to_numpy(),
-        stack,
-        degree=4,
-    )
+    # # generate multichannel projection from label centroids
+    # dff = df[df['keeper']]
+    # interm['projected'] = project_stack_from_focal_points(
+    #     dff['centroid-0'].to_numpy(),
+    #     dff['centroid-1'].to_numpy(),
+    #     dff['zi'].to_numpy(),
+    #     stack,
+    #     degree=4,
+    # )
 
     return {
         'pixel_model_id': px_model.model_id,