#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
cell_map
========
Per-channel cell detection, filtering, atlas alignment, and density mapping.
:class:`CellDetector` is the pipeline worker for the CellMap pipeline, initially
developped to analyze immediate early gene expression data from iDISCO+ cleared tissue [Renier2016]_.
.. image:: ../static/cell_abstract_2016.jpg
:target: https://doi.org/10.1016/j.cell.2020.01.028
:width: 300
.. figure:: ../static/CellMap_pipeline.png
iDISCO+ and ClearMap: A Pipeline for Cell Detection, Registration, and
Mapping in Intact Samples Using Light Sheet Microscopy.
`CellDetector` processes a single fluorescence channel through the following steps:
1. **Cell detection** (:meth:`~CellDetector.run_cell_detection`) —
background subtraction, maxima detection, and shape-based watershed
segmentation via :mod:`~ClearMap.ImageProcessing.Experts.Cells`.
Produces a raw cell table (``cells_raw``).
2. **Filtering** (:meth:`~CellDetector.filter_cells`) —
thresholds on size and intensity; optional cortical-surface crust
removal. Produces a filtered table (``cells_filtered``).
3. **Atlas alignment** (:meth:`~CellDetector.atlas_align`) —
transforms coordinates into atlas space via Elastix and annotates
each cell with brain-region labels, hemisphere, and structure volume.
Produces a per-cell Feather table.
4. **Voxelization** (:meth:`~CellDetector.voxelize`) —
rasterises cell positions (with optional intensity weighting) into a
density volume at atlas resolution for group statistics.
5. **Statistics export** (:meth:`~CellDetector.export_collapsed_stats`) —
aggregates per-cell data into a per-structure CSV with counts, average
size, and volumes for both hemispheres.
Configuration
-------------
Parameters are read from the ``cell_map`` config section:
* ``detection.background_correction.diameter``
* ``detection.maxima_detection.shape`` / ``h_max``
* ``detection.shape_detection.threshold``
* ``cell_filtration.thresholds.size`` / ``intensity``
* ``voxelization.radii``
Typical usage
-------------
Instances are normally created via
:class:`~ClearMap.pipeline_orchestrators.experiment_controller.ExperimentController`::
detector = exp_ctrl.get_worker('cell_map', channel='cfos')
detector.run_cell_detection()
detector.filter_cells()
detector.atlas_align()
detector.voxelize()
detector.export_collapsed_stats()
For scripted use without the GUI::
from ClearMap.pipeline_orchestrators.sample_info_management import build_sample_manager
from ClearMap.pipeline_orchestrators.cell_map import CellDetector
sm = build_sample_manager('/path/to/experiment')
det = CellDetector(sm, sm.cfg_coordinator,
channel='cfos', registration_processor=reg)
det.run_cell_detection()
See also
--------
:mod:`ClearMap.ImageProcessing.Experts.Cells` : Low-level detection algorithms.
:class:`~ClearMap.pipeline_orchestrators.registration_orchestrator.RegistrationProcessor` :
Required for atlas alignment.
:mod:`ClearMap.Analysis.Measurements.Voxelization` : Density-map generation.
"""
from __future__ import annotations
import copy
import re
import platform
import warnings
from concurrent.futures.process import BrokenProcessPool
from typing import Optional, TYPE_CHECKING
import numpy as np
import pandas as pd
# noinspection PyPep8Naming
import ClearMap.IO.IO as clearmap_io
from ClearMap.IO.workspace2 import Workspace2
# noinspection PyPep8Naming
import ClearMap.Alignment.Elastix as elastix
# noinspection PyPep8Naming
import ClearMap.Alignment.Resampling as resampling
# noinspection PyPep8Naming
import ClearMap.ImageProcessing.Experts.Cells as cell_detection
# noinspection PyPep8Naming
import ClearMap.Analysis.Measurements.Voxelization as voxelization
from ClearMap.Utils.exceptions import MissingRequirementException
from ClearMap.Utils.utilities import requires_assets, FilePath, sanitize_n_processes
from ClearMap.config.config_coordinator import ConfigCoordinator
from ClearMap.pipeline_orchestrators.generic_orchestrators import ChannelPipelineOrchestrator
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>, Charly Rousseau <charly.rousseau@icm-institute.org>'
__license__ = 'GPLv3 - GNU General Public License v3 (see LICENSE)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'https://idisco.info'
__download__ = 'https://github.com/ClearAnatomics/ClearMap'
from ClearMap.pipeline_orchestrators.sample_info_management import SampleManager
from ClearMap.pipeline_orchestrators.registration_orchestrator import RegistrationProcessor
if TYPE_CHECKING:
from PyQt5.QtWidgets import QWidget
USE_BINARY_POINTS_FILE = not platform.system().lower().startswith('darwin')
[docs]
class CellDetector(ChannelPipelineOrchestrator):
config_name = 'cell_map'
def __init__(self, sample_manager: SampleManager = None, config_coordinator: ConfigCoordinator = None,
channel=None, registration_processor=None):
super().__init__(config_coordinator)
self.sample_manager: Optional[SampleManager] = None
self.channel: Optional[str] = None
self.registration_processor: Optional[RegistrationProcessor] = None
self.workspace: Optional[Workspace2] = None
self.cell_detection_re = ('Processing block',
re.compile(r'.*?Processing block \d+/\d+.*?\selapsed time:\s\d+:\d+:\d+\.\d+'))
if channel is None:
raise ValueError(f'No channel specified. Please provide a channel name '
f'that matches one in the sample_params file.')
self.setup(sample_manager, channel, registration_processor)
[docs]
def setup(self, sample_manager, channel_name, registration_processor):
self.sample_manager = sample_manager if sample_manager else self.sample_manager
self.channel = channel_name
self.registration_processor = registration_processor
if self.sample_manager.setup_complete:
self.workspace = sample_manager.workspace
self.setup_complete = True
else:
self.setup_complete = False
warnings.warn('SampleManager not setup, CellDetector setup incomplete')
@property
def detected(self):
return self.get('cells', channel=self.channel, asset_sub_type='raw').exists
[docs]
def post_process_cells(self):
self.filter_cells()
self.atlas_align()
self.export_collapsed_stats()
[docs]
def set_voxelization_radii(self, voxelization_radii):
if not isinstance(voxelization_radii, (list, tuple, np.ndarray)):
raise ValueError('voxelization_radius must be a list, tuple or numpy array')
if len(voxelization_radii) != 3:
raise ValueError('voxelization_radius must have three elements (x,y,z)')
if self.config is None:
raise ValueError('CellDetector not properly initialized')
self.patch_channel({'voxelization': {'radii': list(voxelization_radii)}})
[docs]
def list_valid_weighing_columns(self, sub_step: str | None= None):
aligned = sub_step not in ('raw', 'filtered')
cells_df = self.get_coords(coord_type=sub_step, aligned=aligned)
excluded_columns = {'id', 'name', 'order', 'color', 'volume'}
return set(cells_df.columns) - excluded_columns
[docs]
def voxelize(self, sub_step: str | None=None, weights_column=None): # FIXME: add uncrusting ?
"""
Unweighted voxelization (i.e. cell counts)
This will draw a sphere of radius r around each cell and increment the voxel values.
Parameters
----------
sub_step: str | None
If specified, will use the coordinates from the specified sub_step (e.g. 'aligned')
weights_column: str | None
If specified, this column in the cells table will be used to add weights to the
voxelization spheres (e.g. for intensity voxelization).
The column must be present in the cells table.
Returns
-------
coordinates, counts_file_path: np.array, str
"""
if weights_column is not None and weights_column not in self.list_valid_weighing_columns(sub_step=sub_step):
raise ValueError(f'Column {weights_column} is invalid. '
f'Valid options are {self.list_valid_weighing_columns(sub_step)}')
coordinates, cells, voxelization_parameter = self.get_voxelization_params(sub_step=sub_step)
title = 'Voxelisation'
suffix = 'counts'
weights = None
if weights_column:
suffix += f'_{weights_column}'
title += f' weighted by {weights_column}'
weights = self.get_cells_df()[weights_column]
counts_asset = self.get('density', channel=self.channel, asset_sub_type=suffix)
counts_asset.delete(missing_ok=True) # Remove previous counts file if exists
self.set_watcher_step(title)
voxelization.voxelize(coordinates, sink=counts_asset.path, weights=weights, **voxelization_parameter) # WARNING: prange
self.update_watcher_main_progress()
return coordinates, counts_asset.path
@requires_assets([FilePath('density', asset_sub_type='counts')])
def plot_voxelized_counts(self, arrange=True, parent=None):
import ClearMap.Visualization.Plot3d as plot_3d
scale = self.channel_cfg_view('registration')['resampled_resolution']
return plot_3d.plot(self.get_path('density', channel=self.channel, asset_sub_type='counts'),
scale=scale, title='Cell density (voxelized)', lut='flame',
arrange=arrange, parent=parent)
[docs]
def create_test_dataset(self, slicing, debug='debug'):
previous_status = self.workspace.debug
self.workspace.debug = debug
dbg_asset = self.get('stitched', channel=self.channel)
dbg_asset.delete(missing_ok=True) # Remove previous debug asset if exists
self.workspace.debug = previous_status
dbg_path = self.workspace.create_debug('stitched', channel=self.channel, slicing=slicing, debug=debug)
self.update_watcher_main_progress()
return dbg_path
[docs]
def get_voxelization_params(self, sub_step=''):
voxelization_parameter = {
'radius': self.config['voxelization']['radii'],
'verbose': True
}
if self.workspace.debug: # Path will use debug
voxelization_parameter['shape'] = self.get('cells', channel=self.channel, asset_sub_type='shape').shape()
elif self.registration_processor.was_registered:
voxelization_parameter['shape'] = self.get('atlas', channel=self.channel, asset_sub_type='annotation').shape()
else:
voxelization_parameter['shape'] = self.sample_manager.resampled_shape(self.channel)
if sub_step: # Hack to compensate for the fact that the realigned makes no sense in
cells, coordinates = self.get_coords(coord_type=sub_step, aligned=False)
else:
cells, coordinates = self.get_coords(coord_type=None, aligned=True)
return coordinates, cells, voxelization_parameter
[docs]
def get_coords(self, coord_type='filtered', aligned=False):
if coord_type not in ('filtered', 'raw', None):
raise ValueError(f'Coordinate type "{coord_type}" not recognised')
kwargs = {'asset_type': 'cells', 'channel': self.channel, 'asset_sub_type': coord_type}
table_path = self.get_path(**kwargs, extension='.feather')
if not table_path.exists:
table_path = self.get_path(**kwargs)
loaders = {'.feather': pd.read_feather, '.npy': np.load}
table = loaders[table_path.suffix](table_path)
axes = ['xt', 'yt', 'zt'] if aligned else ['x', 'y', 'z']
coordinates = np.array([table[axis] for axis in axes]).T # .T for (n, axes)
return table, coordinates
[docs]
def atlas_align(self):
"""Atlas alignment and annotation """
table, coordinates = self.get_coords(coord_type='filtered')
df = pd.DataFrame({'x': coordinates[:, 0], 'y': coordinates[:, 1], 'z': coordinates[:, 2]})
df['size'] = table['size']
df['source'] = table['source']
if self.registration_processor.was_registered: # FIXME: check if should be registered and raise error
coordinates_transformed = self.transform_coordinates(coordinates)
annotator = self.registration_processor.annotators[self.channel]
ref_channel_cfg = self.get_alignment_ref_channel_reg_cfg()
atlas_resolution = ref_channel_cfg['resampled_resolution']
extra_columns = annotator.get_columns(coordinates_transformed, atlas_resolution) # OPTIMISE: parallel
df = pd.concat([df, extra_columns], axis=1)
else:
warnings.warn('Atlas alignment requires a registered sample.'
'Skipping alignment and using resampled coordinates.')
df.to_feather(self.get_path('cells', channel=self.channel, extension='.feather'))
[docs]
def filter_cells(self, distance_from_surface_px: int = 0):
thresholds = {
'source': self.config['cell_filtration']['thresholds']['intensity'],
'size': self.config['cell_filtration']['thresholds']['size']
}
src_path = self.get_path('cells', channel=self.channel, asset_sub_type='raw')
if not src_path.exists:
raise MissingRequirementException(f'Cell detection not run yet (no file found at "{src_path}"),'
f' cannot filter cells. Please run cell detection first.')
dest_path = self.get_path('cells', channel=self.channel, asset_sub_type='filtered')
cell_detection.filter_cells(source=src_path, sink=dest_path, thresholds=thresholds)
if distance_from_surface_px > 0:
table, filtered_coords = self.get_coords(coord_type='filtered')
uncrusted_coords, mask = self.remove_crust(coordinates=filtered_coords,
threshold=distance_from_surface_px, return_mask=True)
table = table[mask]
clearmap_io.write(dest_path, table) # Overwrite filtered with uncrusted
[docs]
def run_cell_detection(self, tuning=False, save_maxima=False, save_shape=False, save_as_binary_mask=False):
self.workspace.debug = tuning # TODO: use context manager
cell_detection_param = copy.deepcopy(cell_detection.default_cell_detection_parameter)
cell_detection_param['illumination_correction'] = None # WARNING: illumination or illumination_correction
cell_detection_param['background_correction']['shape'] = self.config['detection']['background_correction']['diameter']
cell_detection_param['maxima_detection']['shape'] = self.config['detection']['maxima_detection']['shape']
cell_detection_param['maxima_detection']['h_max'] = self.config['detection']['maxima_detection']['h_max']
cell_detection_param['intensity_detection']['measure'] = ['source']
cell_detection_param['shape_detection']['threshold'] = self.config['detection']['shape_detection']['threshold']
if tuning:
bkg_asset = self.get('cells', channel=self.channel, asset_sub_type='bkg')
bkg_asset.delete(missing_ok=True) # Remove previous background file if exists
cell_detection_param['background_correction']['save'] = str(bkg_asset.path)
shape_asset = self.get('cells', channel=self.channel, asset_sub_type='shape')
if save_shape or tuning:
shape_asset.delete(missing_ok=True)
cell_detection_param['shape_detection']['save'] = str(shape_asset.path)
# if save_as_binary_mask:
# cell_detection_param['shape_detection']['save_dtype'] = 'bool'
if save_maxima:
maxima_asset = self.get('cells', channel=self.channel, asset_sub_type='maxima')
maxima_asset.delete(missing_ok=True)
cell_detection_param['maxima_detection']['save'] = str(maxima_asset.path)
raw_cfg = self.cfg_coordinator.get_config_view(self.config_name)
block_params = raw_cfg['performance']['channels'][self.channel]['detection']['block_processing']
processing_parameter = copy.deepcopy(cell_detection.default_cell_detection_processing_parameter)
# TODO: store as other dict and run .update(**self.extra_detection_params)
processing_parameter.update(processes=sanitize_n_processes(block_params['n_processes']),
size_min=block_params['size_min'], size_max=block_params['size_max'],
overlap=block_params['overlap'], verbose=True)
# TODO: round to processors
n_steps = self.get_n_blocks(self.get('stitched', channel=self.channel).shape()[2])
self.prepare_watcher_for_substep(n_steps, self.cell_detection_re, 'Detecting cells')
try:
dest_asset = self.get('cells', channel=self.channel, asset_sub_type='raw')
dest_asset.delete(missing_ok=True)
cell_detection.detect_cells(self.get_path('stitched', channel=self.channel), dest_asset.path,
cell_detection_parameter=cell_detection_param,
processing_parameter=processing_parameter,
workspace=self.workspace) # WARNING: prange inside multiprocess
# (including array processing and devolve points for vox)
if save_shape and save_as_binary_mask and shape_asset.dtype() != 'bool':
shape_asset.write(np.array(shape_asset.read().astype('bool')))
except BrokenProcessPool as err:
print(f'Cell detection canceled, see: {err}')
return
finally:
self.workspace.debug = False
self.update_watcher_main_progress()
[docs]
def export_as_csv(self):
"""
Export the cell coordinates to csv
.. deprecated:: 2.1
Use :func:`atlas_align` and `export_collapsed_stats` instead.
"""
warnings.warn("export_as_csv is deprecated and will be removed in future versions;"
"please use the new formats from atlas_align and export_collapsed_stats",
DeprecationWarning, 2)
csv_file_path = self.get_path('cells', channel=self.channel, extension='.csv')
self.get_cells_df().to_csv(csv_file_path)
[docs]
def export_collapsed_stats(self, all_regions=True):
df = self.get_cells_df()
collapsed = pd.DataFrame()
relevant_columns = ['id', 'order', 'name', 'hemisphere', 'volume', 'size']
for i in (0, 255): # Split by hemisphere to group by structure and reconcatenate hemispheres after
grouped = df[df['hemisphere'] == i][relevant_columns].groupby(['id'], as_index=False)
tmp = pd.DataFrame()
first = grouped.first()
tmp['Structure ID'] = first['id']
tmp['Structure order'] = first['order']
tmp['Structure name'] = first['name']
tmp['Hemisphere'] = first['hemisphere']
tmp['Structure volume'] = first['volume']
tmp['Cell counts'] = grouped.count()['name']
tmp['Average cell size'] = grouped['size'].mean()['size']
collapsed = pd.concat((collapsed, tmp))
annotator = self.registration_processor.annotators[self.channel]
if all_regions: # Add regions even if they are empty
uniq_ids = np.unique(annotator.atlas)
tmp = pd.DataFrame({'Structure ID': uniq_ids, 'mock': ''})
tmp['Structure name'] = annotator.convert_label(uniq_ids, key='id', value='name')
df_mock = pd.DataFrame({'Hemisphere': [0, 255], 'mock': ''})
tmp = tmp.merge(df_mock, on='mock').drop(columns='mock')
vol_map = annotator.get_lateralised_volume_map(
self.get_alignment_ref_channel_reg_cfg()['resampled_resolution'],
self.get_path('atlas', channel=self.channel, asset_sub_type='hemispheres')
)
tmp['Structure volume'] = tmp.set_index(['Structure ID', 'Hemisphere']).index.map(vol_map.get)
order_map = {int(id_): annotator.find(id_, key='id')['order'] for id_ in uniq_ids}
tmp['Structure order'] = tmp['Structure ID'].map(order_map)
collapsed = tmp.merge(collapsed[['Structure ID', 'Hemisphere', 'Cell counts', 'Average cell size']],
how='left', on=['Structure ID', 'Hemisphere'])
collapsed = collapsed.sort_values(by='Structure ID')
csv_file_path = self.get_path('cells', channel=self.channel, asset_sub_type='stats', extension='.csv')
collapsed.to_csv(csv_file_path, index=False)
[docs]
def plot_cells_3d_scatter_w_atlas_colors(self, raw=False, parent=None):
import ClearMap.Visualization.Qt.Plot3d as qplot_3d
from ClearMap.Visualization.Qt.widgets import Scatter3D
import pyqtgraph as pg
from matplotlib.colors import to_hex
asset_properties = {'channel': self.channel}
if raw:
asset_properties['asset_type'] = 'stitched'
else:
if self.registration_processor.was_registered:
asset_properties['asset_type'] = 'atlas'
asset_properties['asset_sub_type'] = 'reference'
else:
warnings.warn('Dataset not registered, cannot be plotted onto reference atlas. ')
return
asset = self.get(**asset_properties)
requirements = [asset, self.df_path]
if any(not req.exists for req in requirements):
raise MissingRequirementException(f'Cannot plot 3D scatter with atlas colors',
missing_items=[e for e in requirements if not e.exists],
found_items=[e for e in requirements if e.exists])
dv = qplot_3d.plot(asset.path, title=f'{asset_properties["asset_type"].title()} and cells', # FIXME: correct scaling for anisotropic if raw
arrange=False, lut='white', parent=parent)[0]
scatter = pg.ScatterPlotItem()
dv.view.addItem(scatter)
dv.scatter = scatter
df = self.get_cells_df()
if raw:
coordinates = df[['x', 'y', 'z']].values.astype(int)
# coordinates = coordinates * np.array(self.sample_manager.config['resolutions']['raw'])
# coordinates = coordinates.astype(int) # required to match integer z # FIXME: correct scaling for anisotropic
else:
coordinates = df[['xt', 'yt', 'zt']].values.astype(int) # required to match integer z
dv.atlas = self.get('atlas', channel=self.channel, asset_sub_type='annotation').read()
dv.structure_names = self.registration_processor.annotators[self.channel].get_names_map()
if 'hemisphere' in df.columns:
hemispheres = df['hemisphere']
else:
hemispheres = None
if 'id' in df.columns:
unique_ids = np.sort(np.unique(df['id']))
annotator = self.registration_processor.annotators[self.channel]
# color_map = {id_: annotator.find(id_, key='id')['rgb'] for id_ in unique_ids}
color_map = {id_: annotator.convert_label(id_, key='id', value='color_hex_triplet') for id_ in unique_ids}
color_map[0] = np.array(to_hex((1, 0, 0))) # default to red
df['color'] = df['id'].map(color_map)
else:
df['color'] = to_hex((1, 0, 0))
particle_size = self.config['detection']['background_correction']['diameter'][0]
dv.scatter_coords = Scatter3D(coordinates, colors=df['color'].to_list(),
hemispheres=hemispheres,
z_radius=self.machine_config['particle_plot_z_sphere_radius'],
marker_size=max(3, particle_size // 2))
dv.refresh()
return [dv]
@property
def df_path(self):
feather_path = self.get_path('cells', channel=self.channel, extension='.feather')
if feather_path.exists:
return feather_path
else:
return self.get_path('cells', channel=self.channel)
[docs]
def get_cells_df(self):
df_path = self.df_path
if df_path.suffix == '.feather':
return pd.read_feather(df_path)
else:
return pd.DataFrame(np.load(df_path))
@requires_assets([FilePath('cells', asset_sub_type='filtered'), FilePath('stitched')])
def plot_filtered_cells(self, parent=None, smarties=False):
import ClearMap.Visualization.Qt.Plot3d as qplot_3d
from ClearMap.Visualization.Qt.widgets import Scatter3D
import pyqtgraph as pg
_, coordinates = self.get_coords('filtered')
stitched_path = self.get_path('stitched', channel=self.channel)
dv = qplot_3d.plot(stitched_path, title='Stitched and filtered cells', arrange=False,
lut='white', parent=parent)[0]
scatter = pg.ScatterPlotItem()
dv.view.addItem(scatter)
dv.scatter = scatter
dv.scatter_coords = Scatter3D(coordinates, smarties=smarties,
z_radius=self.machine_config['particle_plot_z_sphere_radius'])
dv.refresh()
return [dv]
[docs]
def plot_background_subtracted_img(self):
import ClearMap.Visualization.Plot3d as plot_3d
src = self.get('cells', channel=self.channel, asset_sub_type='raw').open_ro()
coordinates = np.hstack([src[c][:, None] for c in 'xyz'])
p = plot_3d.list_plot_3d(coordinates)
return plot_3d.plot_3d(self.get_path('stitched', channel=self.channel),
view=p, cmap=plot_3d.grays_alpha(alpha=1))
[docs]
def remove_crust(self, coordinates=None, threshold=3, return_mask=False): # TODO: Add inplace option
if coordinates is None:
coordinates = self.get_coords('filtered', aligned=True)
distance_to_surface = self.get('atlas', channel=self.channel, asset_sub_type='distance').read()
# Convert coordinates to integer and insure they are within the distance_to_surface array bounds
int_coordinates = np.floor(coordinates).astype(int) # TODO: check if floor required
xs, ys, zs = int_coordinates.T
xmax, ymax, zmax = distance_to_surface.shape
within_atlas = (xs >= 0) & (xs < xmax) & (ys >= 0) & (ys < ymax) & (zs >= 0) & (zs < zmax)
# Get the distance_to_surface values at the valid coordinates
valid_coordinates = int_coordinates[within_atlas]
dist_values = distance_to_surface[tuple(valid_coordinates.T)]
uncrusted_coordinates = valid_coordinates[dist_values > threshold]
if return_mask:
return uncrusted_coordinates, distance_to_surface[tuple(int_coordinates.T)] > threshold
else:
return uncrusted_coordinates
[docs]
def preview_cell_detection(self, parent: Optional['QWidget'] = None, arrange: bool = True, sync: bool = True) -> list:
import ClearMap.Visualization.Plot3d as plot_3d
sources = [
self.get_path('stitched', channel=self.channel),
self.get_path('cells', channel=self.channel, asset_sub_type='bkg'),
self.get_path('cells', channel=self.channel, asset_sub_type='shape')
]
sources = [s for s in sources if s.exists] # Remove missing files (if not tuning)
if not sources:
raise MissingRequirementException('No files found for preview')
titles = [s.name for s in sources]
luts = ['white', 'white', 'random']
return plot_3d.plot(sources, title=titles, arrange=arrange, sync=sync, lut=luts, parent=parent)
[docs]
def get_n_detected_cells(self):
if self.get('cells', channel=self.channel, asset_sub_type='raw').exists:
_, coords = self.get_coords(coord_type='raw')
return np.max(coords.shape) # TODO: check dimension instead
else:
return 0
[docs]
def get_n_filtered_cells(self):
if self.get('cells', channel=self.channel, asset_sub_type='filtered').exists:
_, coords = self.get_coords(coord_type='filtered')
return np.max(coords.shape) # TODO: check dimension instead
else:
return 0
[docs]
def plot_voxelized_intensities(self, arrange=True):
import ClearMap.Visualization.Plot3d as plot_3d
density_path = self.get_path('density', channel=self.channel, asset_sub_type='intensities')
return plot_3d.plot(density_path, arrange=arrange)
[docs]
def get_n_blocks(self, dim_size):
raw_cfg = self.cfg_coordinator.get_config_view(self.config_name)
perf_params = raw_cfg['performance']['channels'][self.channel]['detection']['block_processing']
blk_size = perf_params['size_max']
overlap = perf_params['overlap']
n_blocks = int(np.ceil((dim_size - blk_size) / (blk_size - overlap) + 1))
return n_blocks
[docs]
def export_to_clearmap1_fmt(self):
"""
ClearMap 1.0 export (will generate the files cells_ClearMap1_intensities, cells_ClearMap1_points_transformed,
cells_ClearMap1_points necessaries to use the analysis script of ClearMap1.
In ClearMap2 the 'cells' file contains already all this information)
In order to align the coordinates when we have right and left hemispheres,
if the orientation of the brain is left, will calculate the new coordinates for the Y axes,
this change will not affect the orientation of the heatmaps, since these are generated from
the ClearMap2 file 'cells'
.. deprecated:: 2.1
Use :func:`atlas_align` and `export_collapsed_stats` instead.
"""
warnings.warn('Method "export_to_clearmap1_fmt" is deprecated and will be removed in future versions;'
'please use the new formats from atlas_align and export_collapsed_stats', DeprecationWarning, 2)
source = self.get('cells', channel=self.channel).read()
clearmap1_format = {'points': ['x', 'y', 'z'],
'points_transformed': ['xt', 'yt', 'zt'],
'intensities': ['source', 'dog', 'background', 'size']}
for sub_type, names in clearmap1_format.items():
sink = self.get_path('cells', channel=self.channel, asset_sub_type=f'ClearMap1{sub_type}')
data = np.array(
[source[name] if name in source.dtype.names else np.full(source.shape[0], np.nan) for name in names]
)
data = data.T
clearmap_io.write(sink, data)
[docs]
def convert_cm2_to_cm2_1_fmt(self):
"""Atlas alignment and annotation """
cells = self.get('cells', channel=self.channel).read()
df = pd.DataFrame({ax: cells[ax] for ax in 'xyz'})
df['size'] = cells['size']
df['source'] = cells['source']
for ax in 'xyz':
df[f'{ax}t'] = cells[f'{ax}t']
df['order'] = cells['order']
df['name'] = cells['name']
coordinates_transformed = np.vstack([cells[f'{ax}t'] for ax in 'xyz']).T
annotator = self.registration_processor.annotators[self.channel]
hemisphere_label = annotator.label_points_hemispheres(coordinates_transformed)
unique_labels = np.sort(df['order'].unique()) # WARNING: ClearMap2 used order as key (now deprecated)
color_map = {lbl: annotator.find(lbl, key='order')['rgb'] for lbl in
unique_labels} # WARNING RGB upper case should give integer but does not work
id_map = {lbl: annotator.find(lbl, key='order')['id'] for lbl in unique_labels}
atlas = self.get('atlas', channel=self.channel, asset_sub_type='annotation').read()
atlas_scale = np.prod(self.get_alignment_ref_channel_reg_cfg['resampled_resolution'])
volumes = {_id: (atlas == _id).sum() * atlas_scale for _id in
id_map.values()} # Volumes need a lookup on ID since the atlas is in ID space
df['id'] = df['order'].map(id_map)
df['hemisphere'] = hemisphere_label
df['color'] = df['order'].map(color_map)
df['volume'] = df['id'].map(volumes)
df.to_feather(self.get_path('cells', channel=self.channel, extension='.feather'))
[docs]
def get_registration_sequence_channels(self, stop_channel='atlas'):
return (self.registration_processor.
get_registration_sequence_channels(self.channel, stop_channel))
