<jupyter_start><jupyter_text>analysis of sequential RNA-FISH for CTP12by <NAME>2021.11.30analysis for dataset:\\10.245.74.158\Chromatin_NAS_1\20211129-P_brain_CTP12-RNA  Load basic functions<jupyter_code>%run "..\..\Startup_py3.py"
sys.path.append(r"..\..\..\..\Documents")

import ImageAnalysis3 as ia
%matplotlib notebook

from ImageAnalysis3 import *
print(os.getpid())

import h5py
from ImageAnalysis3.classes import _allowed_kwds
import ast<jupyter_output>27440<jupyter_text>0. fov parameters<jupyter_code>reload(ia)
reload(classes)
reload(classes.batch_functions)
reload(classes.field_of_view)
reload(io_tools.load)
reload(get_img_info)
reload(visual_tools)
reload(ia.correction_tools)
reload(ia.correction_tools.alignment)
reload(ia.spot_tools.matching)
reload(ia.segmentation_tools.chromosome)
reload(ia.spot_tools.fitting)

fov_param = {'data_folder':r'\\10.245.74.158\Chromatin_NAS_1\20211206-P_brain_CTP12-RNA',
             'save_folder':r'D:\Pu_Temp\20211206-P_brain_CTP12-RNA',
             'experiment_type': 'RNA',
             'num_threads': 16,
             'correction_folder':r'\\10.245.74.158\Chromatin_NAS_0\Corrections\20210621-Corrections_lumencor_from_60_to_50',
             'shared_parameters':{
                     'single_im_size':[50,2048,2048],
                     'distance_zxy': [250, 108, 108],
                     'corr_channels':['750','647'],
                     'num_empty_frames': 0, 
                     'num_buffer_frames':0,
                     'corr_hot_pixel':True,
                     'corr_Z_shift':True,
                     'corr_bleed':True,
                     'min_num_seeds':500,
                     'max_num_seeds': None,
                     'spot_seeding_th': 1500,
                     'normalize_intensity_local':False,
                     'normalize_intensity_background':False,
                     'corr_gaussian_highpass':True
                 }, 
             }

_overwrite = False

for _fov_id in np.arange(41,65):
    # create fov class
    fov = classes.field_of_view.Field_of_View(fov_param, _fov_id=_fov_id,
                                              _color_info_kwargs={
                                                  '_color_filename':'Color_Usage',
                                              }, 
                                              _prioritize_saved_attrs=False,
                                              _save_info_to_file=True,
                                              )
    
    fov.shared_parameters['corr_gaussian_highpass'] = True
    _warp = False
    # 2. Process image into candidate spots
    fov.parallel = True
    fov.rna_ref_id = 0

    fov._process_image_to_spots('rna', 
                                _load_common_reference=True, _load_with_multiple=False,
                                _save_images=True,
                                _warp_images=_warp, 
                                _fit_spots=True,
                                _overwrite_drift=_overwrite, _overwrite_image=_overwrite,
                                _overwrite_spot=_overwrite,
                                _verbose=True)
    
    fov.shared_parameters['corr_gaussian_highpass'] = False

    # 2. Process image into candidate spots
    fov.parallel = False
    fov.protein_ref_id = 0

    fov._process_image_to_spots('protein', 
                                _load_common_reference=True, _load_with_multiple=False,
                                _save_images=True,
                                _warp_images=True, 
                                _fit_spots=False,
                                _overwrite_drift=_overwrite, _overwrite_image=_overwrite,
                                _overwrite_spot=_overwrite,
                                _verbose=True)

    dapi_im = fov._load_dapi_image()
reload(ia.classes.preprocess)
def generate_overlap_crops_3d(crop_size:list or np.ndarray, 
                              image_size:list or np.ndarray, 
                              overlap_size:list or np.ndarray):
    """Function to generate overlapping ImageCrop_3d"""
    _ndim = 3
    crop_size = np.array(crop_size).copy()[:_ndim]
    image_size = np.array(image_size).copy()[:_ndim]
    overlap_size = np.array(overlap_size).copy()[:_ndim]
    
    # check inputs
    if (crop_size > image_size).any():
        raise ValueError(f"crop_size:{crop_size} is larger than image_size:{image_size}")
    if (overlap_size > crop_size).any():
        raise ValueError(f"overlap_size:{overlap_size} is larger than crop_size:{crop_size}")
    
    overlap_size[image_size==crop_size] = 0
    
    # determine crop edges
    num_crops = np.ceil(image_size / crop_size).astype(np.int32)
    
    crop_indices = np.indices(num_crops).reshape(_ndim, -1).transpose()
    
    all_crops = np.ndarray(num_crops, dtype=object)
    
    for _inds in crop_indices:
        left_lim = np.max([_inds * crop_size - overlap_size, np.zeros(_ndim)], axis=0)
        right_lim = np.min([(_inds+1) * crop_size, image_size], axis=0)
        
        _crop = ia.classes.preprocess.ImageCrop_3d(np.array([left_lim, right_lim]).transpose())
        all_crops[tuple(_inds)] = _crop
    
    return all_crops
crops = generate_overlap_crops_3d(
    [50,512,512],
    fov.shared_parameters['single_im_size'],
    [20,128,128],
)
np.ravel(crops)[1].relative_overlap(np.ravel(crops)[0]).array
np.ravel(crops)[1].distance_to_edge([1,2,3])
default_cellpose_kwargs = {
    'diameter': 80, 
    'min_size': 200,
    'stitch_threshold': 0.1,
}

def segmentation_3d_cellpose(raw_im, 
                             pixel_sizes,
                             apply_gaussian=False,
                             gaussian_size=1., 
                             apply_highpass_gaussian=False,
                             highpass_gaussian_kwargs={},
                             cellpose_gpu=True,
                             cellpose_kwargs={},
                             make_plot=True,
                             ):
    """Function to segment single image"""
    if apply_gaussian:
        from scipy.ndimage.filters import gaussian_filter
        _im = gaussian_filter(raw_im, gaussian_size)
    else:
        _im = raw_im.copy()
    
    if apply_highpass_gaussian:
        _im = ia.correction_tools.filters.gaussian_high_pass_filter(
            _im, **highpass_gaussian_kwargs,
        )
    # cellpose segmentation
    import torch
    if torch.cuda.is_available():
        _use_gpu=True
    else:
        _use_gpu=False
    blank_im = np.zeros(np.shape(_im), dtype=np.uint16)
    model = models.Cellpose(gpu=_use_gpu, model_type='nuclei')
    {_k:_v for _k,_v in default_cellpose_kwargs.items()}.update(cellpose_kwargs)

    masks, flows, styles, diams = model.eval(np.stack([blank_im,blank_im,_im], axis=3), 
                                             channels=[0,0],
                                             anisotropy=pixel_sizes[0]/pixel_sizes[1],
                                             do_3D=True,
                                             **default_cellpose_kwargs,
                                             )
    # visualize
    ################
    
    return masks
    
small_masks = []
for _crop in np.ravel(crops):
    _msk = segmentation_3d_cellpose(dapi_im[_crop.to_slices()], 
                                    fov.shared_parameters['distance_zxy'])
    small_masks.append(_msk)
torch.cuda.empty_cache()
?cellpose.models.Cellpose.eval
rescale
_msk = segmentation_3d_cellpose(dapi_im,#[_crop.to_slices()], 
                                fov.shared_parameters['distance_zxy'],
                                
                               cellpose_kwargs={'diameter': 50, 
                                                'rescale':3, 
                                                'resample':False,}
                               )
fov.shared_parameters['distance_zxy']
reload(ia.classes.preprocess)
torch.cuda.empty_cache()
def merge_splitted_masks(crops,
                         masks,
                         image_size, 
                         min_size=1500,
                         merge_per_th=0.,
                         ):
    
    merged_mask = np.zeros(image_size, dtype=np.uint16)
    processed_crops = []
    processed_masks = []
    # loop through masks
    
    for _i, (_crop, _mask) in enumerate(zip(crops, masks)):
        # update mask to prevent sharing labels
        _updated_mask = _mask.copy() + np.max(merged_mask)
        _updated_mask[_updated_mask==np.max(merged_mask)] = 0
        # case 1: if entirely empty, update given area directly
        _target_msk = merged_mask[_crop.to_slices()]
        if (_target_msk==0).all():
            merged_mask[_crop.to_slices()] = _updated_mask
        # case 2: find overlapping area
        else:
            _overlap_crops = [_crop.overlap(_c2) 
                              for _c2 in processed_crops 
                              if _crop.overlap(_c2) is not None]
            # generate overlap masks 
            _overlap_masks = [merge_overlapping_mask(_updated_mask, _crop,
                                                    _m2, _c2, min_size, merge_per_th)
                              for _c2, _m2 in zip(processed_crops, processed_masks) 
                              if _crop.overlap(_c2) is not None]
            # update all first
            merged_mask[_crop.to_slices()] = _updated_mask
            # re-update overlap masks\\
            print(len(_overlap_crops), len(_overlap_masks))
            for _oc, _om in zip(_overlap_crops, _overlap_masks):
                merged_mask[_oc.to_slices()] = _om
                
        # append this processed crop
        processed_crops.append(_crop)
        processed_masks.append(_updated_mask)
    
    
    return merged_mask
mmsk = merge_splitted_masks(np.ravel(crops), small_masks, fov.shared_parameters['single_im_size'])
def merge_overlapping_mask(mask1, crop1,
                           mask2, crop2,
                           min_size=200,
                           merge_per_th=0.,
                           verbose=False,
                           ):
    from scipy.ndimage.measurements import center_of_mass
    if verbose:
        print(f"-- merge overlapping mask")
    
    _o1_crop = crop1.relative_overlap(crop2)
    _o2_crop = crop2.relative_overlap(crop1)
        
    _o1_mask = mask1[_o1_crop.to_slices()]
    _o2_mask = mask2[_o2_crop.to_slices()]
    
    _overlap_mask = np.zeros(np.shape(_o1_mask), dtype=np.uint16)
    
    for _m1 in np.unique(_o1_mask):
        if _m1 > 0:
            _overlap_markers = np.unique(_o2_mask[_o1_mask==_m1])
            # remove small markers
            _o_sizes = np.array([np.sum(_o2_mask==_om) for _om in _overlap_markers])
            _overlap_markers = _overlap_markers[(_overlap_markers>0) * (_o_sizes > min_size)]
            if len(_overlap_markers) == 0:
                if verbose:
                    print(f'--- save uinque mask {_mk}')
                _overlap_mask[_o1_mask==_m1] = _m1
            else:
                # keep the largest sized marker first
                #print(_m1, _overlap_markers)
                #_m1_size = np.sum(mask1==_m1)
                #_m2_sizes = [np.sum(mask1==_m2) for _m2 in _overlap_markers]
                
                _mks = [_m1]+list(_overlap_markers)
                _crop_ref_dicts = {_m1:_o1_mask}
                _crop_ref_dicts.update({_m2:_o2_mask for _m2 in _overlap_markers})
                _ref_dicts = {_m1:[crop1, mask1]}
                _ref_dicts.update({_m2:[crop2, mask2] for _m2 in _overlap_markers})
                # size 
                _sizes = [np.sum(_crop_ref_dicts[_mk]==_mk) for _mk in _mks] #[_m1_size]+list(_m2_sizes)
                # dist
                #_m1_dist = crop1.distance_to_edge(center_of_mass(mask1, labels=mask1, index=[_m1])[0],)
                #_m2_cts = center_of_mass(mask2, labels=mask2, index=_overlap_markers)
                #_m2_dists = [crop2.distance_to_edge(_ct) for _ct in _m2_cts]
                #_dists = [_m1_dist] + _m2_dists
                _dists = [np.sum(mask1==_m1)] + [np.sum(mask2==_m2) for _m2 in _overlap_markers]
                
                # initialize target mask
                _target_mk = _mks[np.argmax(_sizes)]
                _target_dist = _dists[_mks.index(_target_mk)]
                
                #_dists = [_c.distance_to_edge(center_of_mass(_m, labels=_m, index=[_mk])[0],
                #                            ) for _mk, (_c,_m) in _ref_dicts.items()]
                #print('dists', _dists)
                for _i in np.argsort(_sizes)[::-1]:
                    
                    _mk = _mks[_i]
                    _dist = _dists[_i]

                    _ref_mask = _crop_ref_dicts[_mk]
                    _share_mask = (_ref_mask==_mk) * (_overlap_mask==_target_mk)
                    _valid_mask = (_ref_mask==_mk) * (_overlap_mask==0)
                    #print(_mk, np.sum(_share_mask),np.sum(_valid_mask),
                    #      np.sum(_share_mask) / np.sum(_ref_mask==_mk), 
                    #      np.sum(_share_mask) / np.sum(_overlap_mask==_target_mk), )

                    
                    if (np.sum(_share_mask) / np.sum(_ref_mask==_mk) > merge_per_th \
                        or np.sum(_share_mask) / np.sum(_overlap_mask==_target_mk) > merge_per_th):

                        # choose which mask to be used
                        #_target_crop, _target_mask = _ref_dicts[_target_mk]
                        
                        #_target_dist = _target_crop.distance_to_edge(
                        #                    center_of_mass(_target_mask, labels=_target_mask, index=[_target_mk])[0],
                        #                    )
                        #_new_crop, _new_mask = _ref_dicts[_mk]
                        #_new_dist = _new_crop.distance_to_edge(
                        #                    center_of_mass(_new_mask, labels=_new_mask, index=[_mk])[0],
                        #                    )
                        #print(_target_dist, _new_dist)
                        
                        if _target_dist > _dist:
                            if verbose:
                                print(f'--- merge mask {_mk} into {_target_mk}')
                            _overlap_mask[_valid_mask] = _target_mk
                        else:
                            if verbose:
                                print(f'--- merge mask {_target_mk} into {_mk}')
                            _overlap_mask[_valid_mask] = _mk
                            _overlap_mask[_overlap_mask==_target_mk] = _mk
                            # update target mask
                            _target_mk = _mk
                            _target_dist = _dists[_mks.index(_target_mk)]
                    
                    elif np.sum(_valid_mask) > min_size:
                        if verbose:
                            print(f'--- save mask {_mk}')
                        _overlap_mask[_valid_mask] = _mk
                        # update target mask
                        _target_mk = _mk
                        _target_dist = _dists[_mks.index(_target_mk)]
                    else:
                        if verbose:
                            print(f"--- skip mask: {_mk}")
                            
    if verbose:
        print(f"-- {np.sum(np.unique(_overlap_mask) >0)} labels saved.")
    return _overlap_mask#, _o1_mask, _o2_mask
                    
                
                
%%time
ids = [5,6]
_upmsk = small_masks[ids[1]].copy() + small_masks[ids[0]].max()
_upmsk[_upmsk==small_masks[ids[0]].max()] = 0
omsks = merge_overlapping_mask(small_masks[ids[0]], np.ravel(crops)[ids[0]],
                       _upmsk, np.ravel(crops)[ids[1]], 
                               min_size=1500, verbose=True)<jupyter_output>-- merge overlapping mask
{'fraction_non_alphanumeric': 0.096138147566719, 'fraction_numerical': 0.019717425431711146, 'mean_word_length': 1.5522435897435898, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 43, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '610027', 'license': 'MIT', 'n_tokens_mistral': 5153, 'n_tokens_neox': 4832, 'n_words': 725, 'path': 'Tissue_DNA-FISH/CTP12_marker-gene/20211208_mouse-P-brain_CTP12-RNA_franklin.ipynb', 'repo_name': 'shiwei23/Chromatin_Analysis_Scripts', 'size': '15927'}