Merge pull request #17 from templiert/multisem

review of PR #16

Author: trautmane

src/python/janelia_emrp/msem/field_of_view_layout.py

@@ -233,7 +233,7 @@ def __init__(self,
 
     def row_and_col(self,
                     mfov_number: int,
-                    sfov_index_name: str) -> (int, int):
+                    sfov_index_name: str) -> Tuple[int, int]:
         row_offset, col_offset = self.mfov_number_to_offsets[mfov_number]
         sfov_row, sfov_col = self.sfov_index_name_to_row_col[sfov_index_name]
         return sfov_row + row_offset, sfov_col + col_offset

src/python/janelia_emrp/msem/ingestion_ibeammsem/assembly.py

@@ -11,15 +11,14 @@
 import numpy as np
 from dask import bag
 from distributed import Client
-from matplotlib.transforms import Affine2D
-from skimage.io import imread
-from skimage.transform import EuclideanTransform
-
 from janelia_emrp.msem.ingestion_ibeammsem.constant import FACTOR_THUMBNAIL, N_BEAMS
 from janelia_emrp.msem.ingestion_ibeammsem.path import get_image_paths, get_slab_path
 from janelia_emrp.msem.ingestion_ibeammsem.roi import get_mfovs
 from janelia_emrp.msem.ingestion_ibeammsem.xdim import XDim
 from janelia_emrp.msem.ingestion_ibeammsem.xvar import XVar
+from matplotlib.transforms import Affine2D
+from skimage.io import imread
+from skimage.transform import EuclideanTransform
 
 matplotlib.use("tkagg")
 
@@ -29,6 +28,23 @@
     import xarray as xr
 
 
+def get_max_scans(xlog: xr.Dataset) -> int:
+    """Gets the maximum number of scans.
+
+    The xlog is conservatively over-dimensioned upfront along XDim.SCAN
+        to accommodate all anticipated scans.
+
+    The prediction of the number of scans is made by IBEAM-MSEM operators
+        considering the nominal slab thickness
+        and the material removal thickness at every scan.
+
+    Note that scans with strictly negative labels exist,
+        but they are internals of the IBEAM-MSEM process
+        and must not be ingested.
+    """
+    return 1 + xlog[XDim.SCAN].max().values.item()
+
+
 def get_slab_rotation(xlog: xr.Dataset, scan: int, slab: int) -> float:
     """Returns the rotation of a slab, in degrees.
 

src/python/janelia_emrp/msem/ingestion_ibeammsem/id.py

@@ -15,13 +15,22 @@
 
 
 def get_all_magc_ids(xlog: xr.Dataset) -> np.ndarray:
-    """Gets all MagC IDs of the wafer."""
+    """Gets all MagC IDs of the wafer.
+    
+    Note that MagC IDs
+        are not guaranteed to be contiguous, e.g., [0,1,3,5]
+        and do not necessarily start at 0.
+        Therefore use
+            for magc_id in get_all_magc_ids(xlog)
+        instead of
+            for magc_id in range(len(get_all_magc_ids(xlog)))
+    """
     return xlog[XDim.SLAB].values
 
 
 def get_serial_ids(
     xlog: xr.Dataset, magc_ids: list[int] | np.ndarray
-) -> int | None | list[int | None]:
+) -> list[int | None]:
     """Returns the serial IDs of slabs identified by their MagC IDs.
 
     If a magc_id does not have a serial ID, then returns None.

src/python/janelia_emrp/msem/ingestion_ibeammsem/roi.py

@@ -6,7 +6,6 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-
 from janelia_emrp.msem.ingestion_ibeammsem.constant import N_BEAMS
 from janelia_emrp.msem.ingestion_ibeammsem.xdim import XDim
 from janelia_emrp.msem.ingestion_ibeammsem.xvar import XVar
@@ -104,6 +103,20 @@ def get_n_mfovs(xlog: xr.Dataset, scan: int) -> int:
     )
 
 
+def get_max_mfovs_per_slab(xlog: xr.Dataset) -> int:
+    """Gets the maximum number of MFOVs per slab.
+
+    The xlog is dimensioned along XDim.MFOV to fit the slab with the most MFOVs.
+    E.g.: if the largest slab has 24 MFOVs,
+        then the positive labels of the XDim.MFOV are [0,...,23].
+
+    Note that MFOVs with strictly negative IDs exist,
+        but are internals of the IBEAM-MSEM acquisition
+        and are not part of the final dataset.
+    """
+    return 1 + xlog[XDim.MFOV].max().values.item()
+
+
 def get_mfovs(xlog: xr.Dataset, slab: int) -> np.ndarray:
     """Returns the effective MFOV IDs of a slab.
 

src/python/janelia_emrp/msem/msem_to_render.py

@@ -17,11 +17,11 @@
 from janelia_emrp.fibsem.volume_transfer_info import params_to_render_connect
 from janelia_emrp.msem.field_of_view_layout import FieldOfViewLayout, build_mfov_column_group, \
     NINETY_ONE_SFOV_ADJACENT_MFOV_DELTA_Y, NINETY_ONE_SFOV_NAME_TO_ROW_COL
-from janelia_emrp.msem.ingestion_ibeammsem.assembly import get_xys_sfov_and_paths
+from janelia_emrp.msem.ingestion_ibeammsem.assembly import get_xys_sfov_and_paths, get_max_scans
 from janelia_emrp.msem.ingestion_ibeammsem.metrics import get_timestamp
 from janelia_emrp.msem.scan_fit_parameters import ScanFitParameters, \
     build_fit_parameters_path, WAFER_60_61_SCAN_FIT_PARAMETERS
-from janelia_emrp.msem.slab_info import load_slab_info, ContiguousOrderedSlabGroup, MAX_NUMBER_OF_SCANS
+from janelia_emrp.msem.slab_info import load_slab_info, ContiguousOrderedSlabGroup
 from janelia_emrp.root_logger import init_logger
 
 program_name = "msem_to_render.py"
@@ -163,7 +163,7 @@ def import_slab_stacks_for_wafer(render_ws_host: str,
                                  import_magc_slab_list: list[int],
                                  include_scan_list: list[int],
                                  exclude_scan_list: list[int],
-                                 wafer_id: int,
+                                 wafer_id: str,
                                  number_of_slabs_per_render_project: int):
 
     func_name = "import_slab_stacks_for_wafer"
@@ -176,6 +176,9 @@ def import_slab_stacks_for_wafer(render_ws_host: str,
         raise RuntimeError(f"cannot find wafer xlog: {wafer_xlog_path}")
 
     logger.info(f"{func_name}: loading slab info, wafer_id={wafer_id}, number_of_slabs_per_group={number_of_slabs_per_render_project}")
+    
+    n_scans_max = get_max_scans(xlog=xlog)
+    logger.info(f"the maximum number of scans is {n_scans_max}")
 
     slab_group_list = load_slab_info(xlog=xlog,
                                      wafer_id=wafer_id,
@@ -244,7 +247,7 @@ def import_slab_stacks_for_wafer(render_ws_host: str,
                         logger.warning(f'{func_name}: scan {scan} not found for stack {stack}')
             else:
                 # build scan list by looking for first mfov timestamps for all scans and ignoring excluded scans
-                for scan in range(0, MAX_NUMBER_OF_SCANS):
+                for scan in range(0, n_scans_max):
                     first_mfov_scan_timestamp = get_timestamp(xlog=xlog, scan=scan, slab=slab_info.magc_id, mfov=slab_info.first_mfov)
                     if first_mfov_scan_timestamp is not None:
                         if scan not in exclude_scan_list:
@@ -389,7 +392,7 @@ def main(arg_list: List[str]):
     parser.add_argument(
         "--wafer_id",
         help="Wafer identifier (e.g. 60)",
-        type=int,
+        type=str,
         required=True
     )
     parser.add_argument(

src/python/janelia_emrp/msem/slab_info.py

@@ -8,27 +8,31 @@
 from janelia_emrp.msem.ingestion_ibeammsem.id import get_all_magc_ids, get_serial_ids, get_region_ids
 from janelia_emrp.msem.ingestion_ibeammsem.roi import get_mfovs
 
-WAFER_NAME_LEN = 2   # wafers 60 and 61
 SERIAL_NAME_LEN = 3  # 400+ slabs per wafer
 REGION_NAME_LEN = 2  # usually only a few regions per slab, but allow for up to 99
 
-MAX_NUMBER_OF_MFOVS = 100
 MAX_NUMBER_OF_SCANS = 500
 
 @dataclass
 class SlabInfo:
-    wafer_id: int
-    serial_id: int                                     # order in which the slabs were physically cut
+    wafer_id: str
+    serial_id: int
+    """order in which the slabs were physically cut"""
+    magc_id: int = field(compare=False)
+    """magc ID.
+    
+    order in which the slabs were originally defined by the user
+    with the MagFinder plugin in the .magc file
+    """
     region: int
-    magc_id: int = field(compare=False)                # order in which the slabs were originally defined by the user with the MagFinder plugin in the .magc file
     first_mfov: int = field(compare=False)
     last_mfov: int = field(compare=False)
     serial_name: str = field(init=False, repr=False)
     stack_name: str = field(init=False)
 
     def __post_init__(self):
         self.serial_name = f"{self.serial_id:0{SERIAL_NAME_LEN}}"
-        self.stack_name = f"w{self.wafer_id:0{WAFER_NAME_LEN}}_s{self.serial_name}_r{self.region:0{REGION_NAME_LEN}}"
+        self.stack_name = f"w{self.wafer_id}_s{self.serial_name}_r{self.region:0{REGION_NAME_LEN}}"
 
     def build_mfov_position_list(self,
                                  xlog: xarray.Dataset,
@@ -59,16 +63,14 @@ def to_render_project_name(self,
 
 
 def load_slab_info(xlog: xarray.Dataset,
-                   wafer_id: int,
+                   wafer_id: str,
                    number_of_slabs_per_group: int) -> list[ContiguousOrderedSlabGroup]:
 
     magc_ids = get_all_magc_ids(xlog=xlog).tolist()
-    serial_ids = get_serial_ids(xlog=xlog, magc_ids=magc_ids)
 
     slabs: list[SlabInfo] = []
-    for i in range(len(magc_ids)):
-        id_serial=serial_ids[i]
-        slab = magc_ids[i]
+    for slab in magc_ids:
+        id_serial=get_serial_ids(xlog=xlog,magc_ids=[slab])[0]
         mfovs = get_mfovs(xlog=xlog, slab=slab)
         region_ids = get_region_ids(xlog=xlog, slab=slab, mfovs=mfovs)
         id_region = region_ids[0]
@@ -128,7 +130,7 @@ def main(argv: list[str]):
     print(f"loading slab info with wafer_id {argv[2]} and {argv[3]} number_of_slabs_per_group ...")
     number_of_slabs_per_group=int(argv[3])
     slab_groups = load_slab_info(xlog=xlog,
-                                 wafer_id=int(argv[2]),
+                                 wafer_id=argv[2],
                                  number_of_slabs_per_group=number_of_slabs_per_group)
     for slab_group in slab_groups:
         print(f"render project: {slab_group.to_render_project_name(number_of_slabs_per_group)} "