Colocalization cellmap tweak

ClearMap/Analysis/Measurements/maxima_detection.py

@@ -202,3 +202,27 @@ def find_center_of_maxima(source, maxima=None, label=None, verbose=False):
 
     return centers
 
+def label_representatives(A:np.ndarray,eliminate_zero_label=True)->tuple[np.ndarray,...]:
+    """Return the multi-index of one element per non empty labeled region.
+        The 0 labeled region is ignored by default and can be included optionally.
+
+    Parameters
+    ----------
+    A : np.ndarray
+        The labeled regions
+    eliminate_zero_label : bool, optional
+        If the representative for the 0 label is to be eliminated, 
+        by default True
+
+    Returns
+    -------
+    tuple[np.ndarray,...]
+        The multi indices of one representative point per region, the ith 
+        coordinate of the kth point is the kth entry of the ith array in 
+        the tuple.
+    """
+    labels,indices = np.unique(A,return_index=True)
+    if eliminate_zero_label:
+        indices = indices[np.where(labels)]
+    multi_indices = np.unravel_index(indices,A.shape)
+    return multi_indices

ClearMap/Analysis/Measurements/shape_detection.py

@@ -34,57 +34,146 @@
 
 import ClearMap.Utils.Timer as tmr
 import ClearMap.Utils.HierarchicalDict as hdict
-
+from ClearMap.Utils.exceptions import ClearMapValueError
 
 ##############################################################################
 # Cell shape detection
 ##############################################################################
 
-
-def detect_shape(source, seeds, threshold=None, verbose=False, processes=None):
+# version for python >=3.11
+# def labeled_pixels_from_centers(centers,weights,shape):
+#     """label the pixels specify by the centers with the weights.
+
+#     Parameters
+#     ----------
+#     centers : np.ndarray
+#         (n_points,dim) array with points coords, its dtype is the one of weights
+#     weights : np.ndarray
+#         (n_points,) shaped array
+#     shape : tuple
+#         the shape of the output array.
+#     """
+
+#     if centers.shape[0]!=weights.shape[0]:
+#         raise ValueError("The number of weights must equal the number of points.")
+
+#     if len(shape)!=centers.shape[1]:
+#         raise ValueError("Received shape, points with points.shape[1] != len(shape) ")
+
+#     B = np.zeros(shape,dtype = weights.dtype)
+#     transposed = centers.transpose()
+#     B[*(transposed)]=weights
+#     return B
+
+
+def labeled_pixels_from_centers(centers,weights,shape):
+    """label the pixels specify by the centers with the weights.
+
+    This is restricted to 3d shapes for now, a code that will work with python >=3.11
+    is commented out for now
+
+    Parameters
+    ----------
+    centers : np.ndarray
+        (n_points,dim) array with points coords, its dtype is the one of weights
+    weights : np.ndarray
+        (n_points,) shaped array
+    shape : tuple
+        the shape of the output array.
     """
-    Detect object shapes by generating a labeled image from seeds.
-
-    Arguments
-    ---------
+
+    if centers.shape[0]!=weights.shape[0]:
+        raise ValueError("The number of weights must equal the number of points.")
+
+    if len(shape)!=centers.shape[1]:
+        raise ValueError("Received shape, points with points.shape[1] != len(shape) ")
+
+    B = np.zeros(shape,dtype = weights.dtype)
+    xs,ys,zs = centers.transpose()
+    B[xs,ys,zs]=weights
+    return B
+
+def detect_shape(source, seeds, threshold=None, verbose=False, processes=None, as_binary_mask=False, return_sizes=False, seeds_as_labels = False):
+    """Detect object shapes by generating a labeled image from seeds.
+
+    Optionally, the output is replaced by to a mere binary mask and the
+    distinct shapes sizes are also returned.
+
+    Parameters
+    ----------
     source : array, str or Source
         Source image.
     seeds : array, str or Source
-        Cell centers as point coordinates.
+        Cell centers as point coordinates if seeds_as_labels is False. See below otherwise.
     threshold : float or None
         Threshold to determine mask for watershed, pixel below this are
         treated as background. If None, the seeds are expanded indefinitely.
     verbose :bool
         If True, print progress info.
-
+    as_binary_mask : bool, optional
+        If the first output is to be the mask of all shapes, by default False.
+    return_sizes : bool, optional
+        If the sizes of the various shapes are to be returned too, by default False.
+    seeds_as_labels: bool, optional
+        Defaults to Faulse. If True the input seeds is considered to be a labeled image w/
+        the initial basins.
+
     Returns
     -------
     shapes : array
-        Labeled image, where each label indicates an object.
-    """
+        Labeled image, where each label indicates an object. Optionally replaced
+        by shapes>0, see above.
 
+    sizes : (optional) the sizes of the shapes, in the same order are the seeds.
+    """
+    print(f'debug: in detect_shape: threshold={threshold}')
     if verbose:
         timer = tmr.Timer()
         hdict.pprint(head='Shape detection', threshold=threshold)
 
     source = io.as_source(source).array
     seeds = io.as_source(seeds)
-
     mask = None if threshold is None else source > threshold
+    if seeds_as_labels:
+        peaks = seeds
+    else:
+        peaks = labeled_pixels_from_centers(seeds,np.arange(1, seeds.shape[0]+1), source.shape)
+
+    # We check that source has no 0 value otherwise the map source -> -source is not necessarily decreasing, eg for source.dtype=uint16.
+    if np.any(source == 0) and np.issubdtype(source.dtype,np.unsignedinteger):
+        print('For sanity, we shift the source intensity by 1 prior to taking its opposite.')
+        if source.max() < np.ma.minimum_fill_value(source):
+            source+=1
+
+        else:
+            raise ClearMapValueError('An uint array with 0 values will lead to inconsistent results, consider a histogram transform or dtype conversion.')
 
-    peaks = vox.voxelize(seeds, shape=source.shape, weights=np.arange(1, seeds.shape[0]+1), processes=processes).array
     try:
-        shapes = skimage.morphology.watershed(-source, peaks, mask=mask)
+        shapes = skimage.morphology.watershed(-source, peaks, mask=mask, watershed_line=True)
     except AttributeError:
-        shapes = skimage.segmentation.watershed(-source, peaks, mask=mask)
-    # shapes = watershed_ift(-source.astype('uint16'), peaks)
-    # shapes[numpy.logical_not(mask)] = 0
+        shapes = skimage.segmentation.watershed(-source, peaks, mask=mask, watershed_line=True)
+
+    if np.unique(shapes).size != np.unique((peaks if mask is None else peaks*mask)).size:
+        raise RuntimeError(f'watersheding yields unexpected results: the seed number was {np.unique(peaks*mask).size-1}'
+                           + f'and the number of labeled region in output was {np.unique(shapes).size} counting the zero labeled region'
+                           + 'However,' + ( 'there was no zero labeled pixel' if np.count_nonzero(shapes==0)==0  else 'there was some zero labeled pixel') )
 
     if verbose:
         timer.print_elapsed_time('Shape detection')
 
-    return shapes
-
+    if return_sizes:
+        max_label = shapes.max()
+        sizes = find_size(shapes, max_label=max_label)
+
+        if as_binary_mask:
+            return (shapes>0),sizes
+        else:
+            return shapes,sizes
+    else:
+        if as_binary_mask:
+            return shapes>0
+        else:
+            return shapes
 
 def find_size(label, max_label=None, verbose=False):
     """

ClearMap/IO/workspace2.py

@@ -230,15 +230,15 @@ def create_asset(self, type_spec, channel_spec=None, sample_id=None):
         type_spec: TypeSpec
             The type specification of the asset to create.
         channel_spec: ChannelSpec | None
-            The channel to create the asset for.
+            The channel to create the asset fora
         sample_id : str | None
             The sample id to create the asset for. If None, use the workspace sample_id.
             If the workspace sample_id is None, the assets will not use the sample_id
             as prefix in the file name.
 
         Returns
         -------
-
+        asset: Asset
         """
         sample_id = self.sample_id if sample_id is None else sample_id
         asset = Asset(self.directory, type_spec, channel_spec, sample_id=sample_id,
@@ -307,8 +307,9 @@ def get(self, asset_type, channel='current',
 
         if asset_sub_type:
             asset_type += f'_{asset_sub_type}'
-        asset = self.asset_collections[channel][asset_type]
-        if asset is None:
+        if asset_type in self.asset_collections[channel]:
+            asset = self.asset_collections[channel][asset_type]
+        else: # asset is somehow None
             if default == 'closest':
                 warnings.warn('No exact match found. Using partial matching from start.')
                 asset = self.get_closest_matching_asset(asset_type, channel)

ClearMap/ImageProcessing/Experts/Cells.py

@@ -18,11 +18,13 @@
 
 import gc
 import multiprocessing
+import warnings
 
 import numpy as np
 
 import cv2
 import scipy.ndimage.filters as ndf
+import scipy.ndimage as ndi
 
 import ClearMap.IO.IO as clearmap_io
 
@@ -41,6 +43,7 @@
 from ClearMap.ImageProcessing.Experts.utils import initialize_sinks, run_step, print_params
 from ClearMap.ImageProcessing.LocalStatistics import local_percentile
 
+from ClearMap.Utils.exceptions import ClearMapValueError
 ###############################################################################
 # ## Default parameter
 ###############################################################################
@@ -392,11 +395,24 @@ def detect_cells_block(source, parameter=default_cell_detection_parameter, n_thr
         if parameter_maxima['h_max']:  # FIXME: check if source or dog
             centers = md.find_center_of_maxima(source, maxima=maxima, verbose=parameter.get('verbose'))
         else:
-            centers = ap.where(maxima, processes=n_threads).array  # WARNING: prange
+            if parameter_shape:
+                threshold = parameter_shape.get('threshold',0)
+                print(f"masking maxima centers by (dog > {threshold}) for shape detection")
+                mask = dog > threshold
+                maxima = maxima * mask
+
+            # We treat the eventuality of connected components of size>1 in the mask (maxima>0);
+            # the choice of the structure matrix for connectivity could need discussion.
+            # with no further adaptation the maxima_labeling consumes too much memory
+            maxima_labels, _ = ndi.label(maxima, structure=np.ones((3,)*3,dtype='bool'))
+            centers = np.vstack(md.label_representatives(maxima_labels)).transpose()
+            # we could come back to the ancient version
+            # centers = ap.where(maxima, processes=n_threads).array 
         del maxima
 
         # correct for valid region
         if valid:
+            print('Filtering centers for correct block processing.')
             ids = np.ones(len(centers), dtype=bool)
             for c, l, u in zip(centers.T, valid_lower, valid_upper):
                 ids = np.logical_and(ids, np.logical_and(l <= c, c < u))
@@ -406,16 +422,24 @@ def detect_cells_block(source, parameter=default_cell_detection_parameter, n_thr
     else:
         results = ()
 
-    # WARNING: sd.detect_shape uses prange
     # cell shape detection  # FIXME: may use centers without assignment
-    shape = run_step('shape_detection', dog, sd.detect_shape, remove_previous_result=True, **default_step_params,
-                     args=[centers], extra_kwargs={'verbose': parameter.get('verbose'), 'processes': n_threads})
     if parameter_shape:
-        # size detection
-        max_label = centers.shape[0]
-        sizes = sd.find_size(shape, max_label=max_label)
+        try:
+            parser = (lambda t: t[0]>0)
+            shape, sizes = run_step('shape_detection', dog, sd.detect_shape, remove_previous_result=True, **default_step_params,
+            args=[centers], presave_parser=parser, extra_kwargs={'verbose': parameter.get('verbose'), 'processes': n_threads, 'return_sizes': True})
+        except ClearMapValueError as err:
+            if str(err) == 'An uint array with 0 values will lead to inconsistent results, consider a histogram transform or dtype conversion.':
+                warnings.warn('This block is likely to contain corrupted data, an empty output will be provided for this block.')
+                results = (centers[:0],)
+                sizes = np.array([])
+                shape = None
+            else:
+                raise err
+
+
+
         valid = sizes > 0
-
         results += (sizes,)
     else:
         valid = None
@@ -433,11 +457,15 @@ def detect_cells_block(source, parameter=default_cell_detection_parameter, n_thr
 
         for m in measure:
             if shape is not None:
+                max_label = centers.shape[0]
                 intensity = sd.find_intensity(steps_to_measure[m], label=shape,
-                                              max_label=max_label, **parameter_intensity)
+                                            max_label=max_label, **parameter_intensity)
             else:  # WARNING: prange but me.measure_expression not parallel since processes=1
+                # FIXME : How can r be defined in this branch ???
                 intensity = me.measure_expression(steps_to_measure[m], centers, search_radius=r,
-                                                  **parameter_intensity, processes=1, verbose=False)
+
+                                                    **parameter_intensity, processes=1, verbose=False)
+
             results += (intensity,)
 
         if parameter.get('verbose'):

ClearMap/ImageProcessing/Experts/utils.py

@@ -11,7 +11,10 @@ def initialize_sinks(cell_detection_parameter, shape, order):
         if isinstance(par, dict):
             filename = par.get('save')
             if filename:
-                ap.initialize_sink(filename, shape=shape, order=order, dtype='float')
+                dtype = par.get('save_dtype')
+                if dtype is None:
+                    dtype='float'
+                ap.initialize_sink(filename, shape=shape, order=order, dtype=dtype)
 
 
 def print_params(step_params, param_key, prefix, verbose):
@@ -25,25 +28,33 @@ def print_params(step_params, param_key, prefix, verbose):
 
 
 def run_step(param_key, previous_result, step_function, args=(), remove_previous_result=False,
-             extra_kwargs=None, parameter=None, steps_to_measure=None, prefix='',
+             presave_parser=None, extra_kwargs=None, parameter=None, steps_to_measure=None, prefix='',
              base_slicing=None, valid_slicing=None):
     if steps_to_measure is None:
         steps_to_measure = {}
     if parameter is None:
         parameter = {}
     if extra_kwargs is None:
         extra_kwargs = {}
-
     step_param = parameter.get(param_key)
     if step_param:
         step_param, timer = print_params(step_param, param_key, prefix, parameter.get('verbose'))
-
         save = step_param.pop('save', None)  # FIXME: check if always goes before step_function call
+        save_dtype = step_param.pop('save_dtype', None)  
         result = step_function(previous_result, *args, **{**step_param, **extra_kwargs})
 
         if save:
-            save = clearmap_io.as_source(save)
-            save[base_slicing] = result[valid_slicing]
+            if save_dtype is None:
+                save = clearmap_io.as_source(save)
+            else:                    
+                save = clearmap_io.as_source(save,dtype=save_dtype)
+
+            to_save=presave_parser(result)
+            # usefull sanity check
+            if save_dtype=='bool':
+                save[base_slicing] = (to_save[valid_slicing] > 0)
+            else:
+                save[base_slicing] = to_save[valid_slicing]
 
         if parameter.get('verbose'):
             timer.print_elapsed_time(param_key.title())