# HG changeset patch
# User goeckslab
# Date 1663358754 0
# Node ID 96d0d969ebc99e45e4f46e10f748b368c70efd27
# Parent 41e8efe8df43665738c295dde051a5d277cae65a
planemo upload for repository https://github.com/goeckslab/tools-mti/tree/main/tools/s3segmenter commit 0f4f17235c5961c2fd3d4c30180507f66214c11d
diff -r 41e8efe8df43 -r 96d0d969ebc9 S3segmenter.py
--- a/S3segmenter.py Fri Mar 11 23:37:49 2022 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,592 +0,0 @@
-import matplotlib.pyplot as plt
-import tifffile
-import os
-import numpy as np
-from skimage import io as skio
-from scipy.ndimage import *
-import scipy.ndimage as ndi
-from skimage.morphology import *
-from skimage.morphology import extrema
-from skimage import morphology
-from skimage.measure import regionprops
-from skimage.transform import resize
-from skimage.filters import gaussian, threshold_otsu, threshold_local
-from skimage.feature import peak_local_max
-from skimage.color import label2rgb
-from skimage.io import imsave,imread
-from skimage.segmentation import clear_border, watershed, find_boundaries
-from scipy.ndimage.filters import uniform_filter
-from os.path import *
-from os import listdir, makedirs, remove
-import pickle
-import shutil
-import fnmatch
-import cv2
-import sys
-import argparse
-import re
-import copy
-import datetime
-from joblib import Parallel, delayed
-from rowit import WindowView, crop_with_padding_mask
-from save_tifffile_pyramid import save_pyramid
-import subprocess
-import ome_types
-
-
-def imshowpair(A,B):
- plt.imshow(A,cmap='Purples')
- plt.imshow(B,cmap='Greens',alpha=0.5)
- plt.show()
-
-
-def imshow(A):
- plt.imshow(A)
- plt.show()
-
-def overlayOutline(outline,img):
- img2 = img.copy()
- stacked_img = np.stack((img2,)*3, axis=-1)
- stacked_img[outline > 0] = [65535, 0, 0]
- imshowpair(img2,stacked_img)
-
-def normI(I):
- Irs=resize(I,(I.shape[0]//10,I.shape[1]//10) );
- p1 = np.percentile(Irs,10);
- J = I-p1;
- p99 = np.percentile(Irs,99.99);
- J = J/(p99-p1);
- return J
-
-def contour_pm_watershed(
- contour_pm, sigma=2, h=0, tissue_mask=None,
- padding_mask=None, min_area=None, max_area=None
-):
- if tissue_mask is None:
- tissue_mask = np.ones_like(contour_pm)
- padded = None
- if padding_mask is not None and np.any(padding_mask == 0):
- contour_pm, padded = crop_with_padding_mask(
- contour_pm, padding_mask, return_mask=True
- )
- tissue_mask = crop_with_padding_mask(
- tissue_mask, padding_mask
- )
-
- maxima = peak_local_max(
- extrema.h_maxima(
- ndi.gaussian_filter(np.invert(contour_pm), sigma=sigma),
- h=h
- ),
- indices=False,
- footprint=np.ones((3, 3))
- )
- maxima = label(maxima).astype(np.int32)
-
- # Passing mask into the watershed function will exclude seeds outside
- # of the mask, which gives fewer and more accurate segments
- maxima = watershed(
- contour_pm, maxima, watershed_line=True, mask=tissue_mask
- ) > 0
-
- if min_area is not None and max_area is not None:
- maxima = label(maxima, connectivity=1).astype(np.int32)
- areas = np.bincount(maxima.ravel())
- size_passed = np.arange(areas.size)[
- np.logical_and(areas > min_area, areas < max_area)
- ]
- maxima *= np.isin(maxima, size_passed)
- np.greater(maxima, 0, out=maxima)
-
- if padded is None:
- return maxima.astype(np.bool)
- else:
- padded[padded == 1] = maxima.flatten()
- return padded.astype(np.bool)
-
-def S3AreaSegmenter(nucleiPM, images, TMAmask, threshold,fileprefix,outputPath):
- nucleiCenters = nucleiPM[:,:,0]
- TMAmask= (nucleiCenters>np.amax(nucleiCenters)*0.8)*TMAmask
- area = []
- area.append(np.sum(np.sum(TMAmask)))
- for iChan in range(len(images)):
- image_gauss = gaussian(resize(images[iChan,:,:],(int(0.25*images.shape[1]),int(0.25*images.shape[2]))),1)
- if threshold ==-1:
- threshold = threshold_otsu(image_gauss)
- mask=resize(image_gauss>threshold,(images.shape[1],images.shape[2]),order = 0)*TMAmask
- area.append(np.sum(np.sum(mask)))
- os.mk
- np.savetxt(outputPath + os.path.sep + fileprefix + '_area.csv',(np.transpose(np.asarray(area))),fmt='%10.5f')
- return TMAmask
-
-def getMetadata(path,commit):
- with tifffile.TiffFile(path) as tif:
- if not tif.ome_metadata:
- try:
- x_res_tag = tif.pages[0].tags['XResolution'].value
- y_res_tag = tif.pages[0].tags['YResolution'].value
- physical_size_x = x_res_tag[0] / x_res_tag[1]
- physical_size_y = y_res_tag[0] / y_res_tag[1]
- except KeyError:
- physical_size_x = 1
- physical_size_y = 1
- metadata_args = dict(
- pixel_sizes=(physical_size_y, physical_size_x),
- pixel_size_units=('µm', 'µm'),
- software= 's3segmenter v' + commit
- )
- else:
- metadata=ome_types.from_xml(tif.ome_metadata)
- metadata = metadata.images[0].pixels
- metadata_args = dict(
- pixel_sizes=(metadata.physical_size_y,metadata.physical_size_x),
- pixel_size_units=('µm', 'µm'),
- software= 's3segmenter v' + commit
- )
- return metadata_args
-
-def S3NucleiBypass(nucleiPM,nucleiImage,logSigma,TMAmask,nucleiFilter,nucleiRegion):
- foregroundMask = nucleiPM
- P = regionprops(foregroundMask, nucleiImage)
- prop_keys = ['mean_intensity', 'label','area']
- def props_of_keys(prop, keys):
- return [prop[k] for k in keys]
-
- mean_ints, labels, areas = np.array(Parallel(n_jobs=6)(delayed(props_of_keys)(prop, prop_keys)
- for prop in P
- )
- ).T
- del P
- maxArea = (logSigma[1]**2)*3/4
- minArea = (logSigma[0]**2)*3/4
- passed = np.logical_and(areas > minArea, areas < maxArea)
-
- foregroundMask *= np.isin(foregroundMask, labels[passed])
- np.greater(foregroundMask, 0, out=foregroundMask)
- foregroundMask = label(foregroundMask, connectivity=1).astype(np.int32)
- return foregroundMask
-
-def S3NucleiSegmentationWatershed(nucleiPM,nucleiImage,logSigma,TMAmask,nucleiFilter,nucleiRegion):
- nucleiContours = nucleiPM[:,:,1]
- nucleiCenters = nucleiPM[:,:,0]
- mask = resize(TMAmask,(nucleiImage.shape[0],nucleiImage.shape[1]),order = 0)>0
- if nucleiRegion == 'localThreshold' or nucleiRegion == 'localMax':
- Imax = peak_local_max(extrema.h_maxima(255-nucleiContours,logSigma[0]),indices=False)
- Imax = label(Imax).astype(np.int32)
- foregroundMask = watershed(nucleiContours, Imax, watershed_line=True)
- P = regionprops(foregroundMask, np.amax(nucleiCenters) - nucleiCenters - nucleiContours)
- prop_keys = ['mean_intensity', 'label','area']
- def props_of_keys(prop, keys):
- return [prop[k] for k in keys]
-
- mean_ints, labels, areas = np.array(Parallel(n_jobs=6)(delayed(props_of_keys)(prop, prop_keys)
- for prop in P
- )
- ).T
- del P
- maxArea = (logSigma[1]**2)*3/4
- minArea = (logSigma[0]**2)*3/4
- passed = np.logical_and.reduce((
- np.logical_and(areas > minArea, areas < maxArea),
- np.less(mean_ints, 50)
- ))
-
- foregroundMask *= np.isin(foregroundMask, labels[passed])
- np.greater(foregroundMask, 0, out=foregroundMask)
- foregroundMask = label(foregroundMask, connectivity=1).astype(np.int32)
-
- else:
- if len(logSigma)==1:
- nucleiDiameter = [logSigma*0.5, logSigma*1.5]
- else:
- nucleiDiameter = logSigma
- logMask = nucleiCenters > 150
-
- win_view_setting = WindowView(nucleiContours.shape, 2000, 500)
-
- nucleiContours = win_view_setting.window_view_list(nucleiContours)
- padding_mask = win_view_setting.padding_mask()
- mask = win_view_setting.window_view_list(mask)
-
- maxArea = (logSigma[1]**2)*3/4
- minArea = (logSigma[0]**2)*3/4
-
- foregroundMask = np.array(
- Parallel(n_jobs=6)(delayed(contour_pm_watershed)(
- img, sigma=logSigma[1]/30, h=logSigma[1]/30, tissue_mask=tm,
- padding_mask=m, min_area=minArea, max_area=maxArea
- ) for img, tm, m in zip(nucleiContours, mask, padding_mask))
- )
-
- del nucleiContours, mask
-
- foregroundMask = win_view_setting.reconstruct(foregroundMask)
- foregroundMask = label(foregroundMask, connectivity=1).astype(np.int32)
-
- if nucleiFilter == 'IntPM':
- int_img = nucleiCenters
- elif nucleiFilter == 'Int':
- int_img = nucleiImage
-
- print(' ', datetime.datetime.now(), 'regionprops')
- P = regionprops(foregroundMask, int_img)
-
- def props_of_keys(prop, keys):
- return [prop[k] for k in keys]
-
- prop_keys = ['mean_intensity', 'area', 'solidity', 'label']
- mean_ints, areas, solidities, labels = np.array(
- Parallel(n_jobs=6)(delayed(props_of_keys)(prop, prop_keys)
- for prop in P
- )
- ).T
- del P
-
- MITh = threshold_otsu(mean_ints)
-
- minSolidity = 0.8
-
- passed = np.logical_and.reduce((
- np.greater(mean_ints, MITh),
- np.logical_and(areas > minArea, areas < maxArea),
- np.greater(solidities, minSolidity)
- ))
-
- # set failed mask label to zero
- foregroundMask *= np.isin(foregroundMask, labels[passed])
-
- np.greater(foregroundMask, 0, out=foregroundMask)
- foregroundMask = label(foregroundMask, connectivity=1).astype(np.int32)
-
- return foregroundMask
-
-def bwmorph(mask,radius):
- mask = np.array(mask,dtype=np.uint8)
- #labels = label(mask)
- background = nucleiMask == 0
- distances, (i, j) = distance_transform_edt(background, return_indices=True)
- cellMask = nucleiMask.copy()
- finalmask = background & (distances <= radius)
- cellMask[finalmask] = nucleiMask[i[finalmask], j[finalmask]]
-
-# imshowpair(cellMask,mask)
- return cellMask
-# imshow(fg)
-# fg = cv2.dilate(mask,ndimage.generate_binary_structure(2, 2))
-# bg = 1-fg-mask
-# imshowpair(bg,mask)
-
-def S3CytoplasmSegmentation(nucleiMask,cyto,mask,cytoMethod='distanceTransform',radius = 5):
- mask = (nucleiMask + resize(mask,(nucleiMask.shape[0],nucleiMask.shape[1]),order=0))>0
- gdist = distance_transform_edt(1-(nucleiMask>0))
- if cytoMethod == 'distanceTransform':
- mask = np.array(mask,dtype=np.uint32)
- markers= nucleiMask
- elif cytoMethod == 'hybrid':
- cytoBlur = gaussian(cyto,2)
- c1 = uniform_filter(cytoBlur, 3, mode='reflect')
- c2 = uniform_filter(cytoBlur*cytoBlur, 3, mode='reflect')
- grad = np.sqrt(c2 - c1*c1)*np.sqrt(9./8)
- grad[np.isnan(grad)]=0
- gdist= np.sqrt(np.square(grad) + 0.000001*np.amax(grad)/np.amax(gdist)*np.square(gdist))
- bg = binary_erosion(np.invert(mask),morphology.selem.disk(radius, np.uint8))
- markers=nucleiMask.copy()
- markers[bg==1] = np.amax(nucleiMask)+1
- markers = label(markers>0,connectivity=1)
- mask = np.ones(nucleiMask.shape)
- del bg
- elif cytoMethod == 'ring':
-# mask =np.array(bwmorph(nucleiMask,radius)*mask,dtype=np.uint32)>0
- mask =np.array(bwmorph(nucleiMask,radius),dtype=np.uint32)>0
- markers= nucleiMask
-
- cellMask =clear_border(watershed(gdist,markers,watershed_line=True))
- del gdist, markers, cyto
- cellMask = np.array(cellMask*mask,dtype=np.uint32)
-
- finalCellMask = np.zeros(cellMask.shape,dtype=np.uint32)
- P = regionprops(label(cellMask>0,connectivity=1),nucleiMask>0,cache=False)
- count=0
- for props in P:
- if props.max_intensity>0 :
- count += 1
- yi = props.coords[:, 0]
- xi = props.coords[:, 1]
- finalCellMask[yi, xi] = count
- nucleiMask = np.array(nucleiMask>0,dtype=np.uint32)
- nucleiMask = finalCellMask*nucleiMask
- cytoplasmMask = np.subtract(finalCellMask,nucleiMask)
- return cytoplasmMask,nucleiMask,finalCellMask
-
-def exportMasks(mask,image,outputPath,filePrefix,fileName,commit,metadata_args,saveFig=True,saveMasks = True):
- outputPath =outputPath + os.path.sep + filePrefix
- if not os.path.exists(outputPath):
- os.makedirs(outputPath)
- previewPath = outputPath + os.path.sep + 'qc'
- if not os.path.exists(previewPath):
- os.makedirs(previewPath)
-
- if saveMasks ==True:
- save_pyramid(
- mask,
- outputPath + os.path.sep + fileName + '.ome.tif',
- channel_names=fileName,
- is_mask=True,
- **metadata_args
- )
- if saveFig== True:
- mask=np.uint8(mask>0)
- edges = find_boundaries(mask,mode = 'outer')
- stacked_img=np.stack((np.uint16(edges)*np.amax(image),image),axis=0)
- save_pyramid(
- stacked_img,
- previewPath + os.path.sep + fileName + 'Outlines.ome.tif',
- channel_names=[f'{fileName} outlines', 'Segmentation image'],
- is_mask=False,
- **metadata_args
- )
-
-def S3punctaDetection(spotChan,mask,sigma,SD):
- Ilog = -gaussian_laplace(np.float32(spotChan),sigma)
- tissueMask = spotChan >0
- fgm=peak_local_max(Ilog, indices=False,footprint=np.ones((3, 3)))*tissueMask
- test=Ilog[fgm==1]
- med = np.median(test)
- mad = np.median(np.absolute(test - med))
- thresh = med + 1.4826*SD*mad
- return (Ilog>thresh)*fgm*(mask>0)
-
-
-
-if __name__ == '__main__':
- parser=argparse.ArgumentParser()
- parser.add_argument("--imagePath")
- parser.add_argument("--contoursClassProbPath",default ='')
- parser.add_argument("--nucleiClassProbPath",default ='')
- parser.add_argument("--stackProbPath",default ='')
- parser.add_argument("--outputPath")
- parser.add_argument("--dearrayPath")
- parser.add_argument("--maskPath")
- parser.add_argument("--probMapChan",type = int, default = -1)
- parser.add_argument("--mask",choices=['TMA', 'tissue','none'],default = 'tissue')
- parser.add_argument("--crop",choices=['interactiveCrop','autoCrop','noCrop','dearray','plate'], default = 'noCrop')
- parser.add_argument("--cytoMethod",choices=['hybrid','distanceTransform','bwdistanceTransform','ring'],default = 'distanceTransform')
- parser.add_argument("--nucleiFilter",choices=['IntPM','LoG','Int','none'],default = 'IntPM')
- parser.add_argument("--nucleiRegion",choices=['watershedContourDist','watershedContourInt','watershedBWDist','dilation','localThreshold','localMax','bypass','pixellevel'], default = 'watershedContourInt')
- parser.add_argument("--pixelThreshold",type = float, default = -1)
- parser.add_argument("--segmentCytoplasm",choices = ['segmentCytoplasm','ignoreCytoplasm'],default = 'ignoreCytoplasm')
- parser.add_argument("--cytoDilation",type = int, default = 5)
- parser.add_argument("--logSigma",type = int, nargs = '+', default = [3, 60])
- parser.add_argument("--CytoMaskChan",type=int, nargs = '+', default=[2])
- parser.add_argument("--pixelMaskChan",type=int, nargs = '+', default=[2])
- parser.add_argument("--TissueMaskChan",type=int, nargs = '+', default=0)
- parser.add_argument("--detectPuncta",type=int, nargs = '+', default=[0])
- parser.add_argument("--punctaSigma", nargs = '+', type=float, default=[0])
- parser.add_argument("--punctaSD", nargs = '+', type=float, default=[4])
- parser.add_argument("--saveMask",action='store_false')
- parser.add_argument("--saveFig",action='store_false')
- args = parser.parse_args()
-
- imagePath = args.imagePath
- outputPath = args.outputPath
- nucleiClassProbPath = args.nucleiClassProbPath
- contoursClassProbPath = args.contoursClassProbPath
- stackProbPath = args.stackProbPath
- maskPath = args.maskPath
-
- commit = '1.3.11'#subprocess.check_output(['git', 'describe', '--tags']).decode('ascii').strip()
- metadata = getMetadata(imagePath,commit)
-
- fileName = os.path.basename(imagePath)
- filePrefix = fileName[0:fileName.index('.')]
-
- # convert 1-based indexing to 0-based indexing
- CytoMaskChan = args.CytoMaskChan
- CytoMaskChan[:] = [number - 1 for number in CytoMaskChan]
- pixelMaskChan = args.pixelMaskChan
- pixelMaskChan[:] = [number - 1 for number in pixelMaskChan]
-
-
- if not os.path.exists(outputPath):
- os.makedirs(outputPath)
-
-
- # get channel used for nuclei segmentation
-
- if len(contoursClassProbPath)>0:
- legacyMode = 1
- probPrefix = os.path.basename(contoursClassProbPath)
- nucMaskChan = int(probPrefix.split('ContoursPM_')[1].split('.')[0])
- elif len(stackProbPath)>0:
- legacyMode = 0
- probPrefix = os.path.basename(stackProbPath)
- else:
- print('NO PROBABILITY MAP PROVIDED')
-
- if args.probMapChan==-1:
- print('Using first channel by default!')
- nucMaskChan = 0
- else:
- nucMaskChan = args.probMapChan
- nucMaskChan = nucMaskChan -1 #convert 1-based indexing to 0-based indexing
-
- if args.TissueMaskChan==0:
- TissueMaskChan = copy.copy(CytoMaskChan)
- TissueMaskChan.append(nucMaskChan)
- else:
- TissueMaskChan = args.TissueMaskChan[:]
- TissueMaskChan[:] = [number - 1 for number in TissueMaskChan]#convert 1-based indexing to 0-based indexing
- TissueMaskChan.append(nucMaskChan)
-
- #crop images if needed
- if args.crop == 'interactiveCrop':
- nucleiCrop = tifffile.imread(imagePath,key = nucMaskChan)
- r=cv2.selectROI(resize(nucleiCrop,(nucleiCrop.shape[0] // 30, nucleiCrop.shape[1] // 30)))
- cv2.destroyWindow('select')
- rect=np.transpose(r)*30
- PMrect= [rect[1], rect[0], rect[3], rect[2]]
- nucleiCrop = nucleiCrop[int(rect[1]):int(rect[1]+rect[3]), int(rect[0]):int(rect[0]+rect[2])]
- elif args.crop == 'noCrop' or args.crop == 'dearray' or args.crop == 'plate':
- nucleiCrop = tifffile.imread(imagePath,key = nucMaskChan)
- rect = [0, 0, nucleiCrop.shape[0], nucleiCrop.shape[1]]
- PMrect= rect
- elif args.crop == 'autoCrop':
- nucleiCrop = tifffile.imread(imagePath,key = nucMaskChan)
- rect = [np.round(nucleiCrop.shape[0]/3), np.round(nucleiCrop.shape[1]/3),np.round(nucleiCrop.shape[0]/3), np.round(nucleiCrop.shape[1]/3)]
- PMrect= rect
- nucleiCrop = nucleiCrop[int(rect[0]):int(rect[0]+rect[2]), int(rect[1]):int(rect[1]+rect[3])]
-
- if legacyMode ==1:
- nucleiProbMaps = tifffile.imread(nucleiClassProbPath,key=0)
- nucleiPM = nucleiProbMaps[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- nucleiProbMaps = tifffile.imread(contoursClassProbPath,key=0)
- PMSize = nucleiProbMaps.shape
- nucleiPM = np.dstack((nucleiPM,nucleiProbMaps[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]))
- else:
- nucleiProbMaps = imread(stackProbPath)
- if len(nucleiProbMaps.shape)==2:
- nucleiPM = nucleiProbMaps[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- else:
- nucleiPM = nucleiProbMaps[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3]),:]
- PMSize = nucleiProbMaps.shape
-
- # mask the core/tissue
- if args.crop == 'dearray':
- TMAmask = tifffile.imread(maskPath)
- elif args.crop =='plate':
- TMAmask = np.ones(nucleiCrop.shape)
-
- else:
- tissue = np.empty((len(TissueMaskChan),nucleiCrop.shape[0],nucleiCrop.shape[1]),dtype=np.uint16)
- count=0
- if args.crop == 'noCrop':
- for iChan in TissueMaskChan:
- tissueCrop =tifffile.imread(imagePath,key=iChan)
- tissue_gauss = gaussian(tissueCrop,1)
- #tissue_gauss[tissue_gauss==0]=np.nan
- tissue[count,:,:] =np.log2(tissue_gauss+1)>threshold_otsu(np.log2(tissue_gauss+1))
- count+=1
- else:
- for iChan in TissueMaskChan:
- tissueCrop = tifffile.imread(imagePath,key=iChan)
- tissue_gauss = gaussian(tissueCrop[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])],1)
- tissue[count,:,:] = np.log2(tissue_gauss+1)>threshold_otsu(np.log2(tissue_gauss+1))
- count+=1
- TMAmask = np.max(tissue,axis = 0)
-
-
- del tissue_gauss, tissue
-
- # nuclei segmentation
- if args.nucleiRegion == 'pixellevel':
- pixelTissue = np.empty((len(pixelMaskChan),nucleiCrop.shape[0],nucleiCrop.shape[1]),dtype=np.uint16)
- count=0
- for iChan in pixelMaskChan:
- pixelCrop = tifffile.imread(imagePath,key=iChan)
- pixelTissue[count,:,:] = pixelCrop[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- count+=1
- nucleiMask = S3AreaSegmenter(nucleiPM, pixelTissue, TMAmask,args.pixelThreshold,filePrefix,outputPath)
- elif args.nucleiRegion == 'bypass':
- nucleiMask = S3NucleiBypass(nucleiPM,nucleiCrop,args.logSigma,TMAmask,args.nucleiFilter,args.nucleiRegion)
- else:
- nucleiMask = S3NucleiSegmentationWatershed(nucleiPM,nucleiCrop,args.logSigma,TMAmask,args.nucleiFilter,args.nucleiRegion)
- del nucleiPM
- # cytoplasm segmentation
- if args.segmentCytoplasm == 'segmentCytoplasm':
- count =0
- if args.crop == 'noCrop' or args.crop == 'dearray' or args.crop == 'plate':
- cyto=np.empty((len(CytoMaskChan),nucleiCrop.shape[0],nucleiCrop.shape[1]),dtype=np.uint16)
- for iChan in CytoMaskChan:
- cyto[count,:,:] = tifffile.imread(imagePath, key=iChan)
- count+=1
- elif args.crop =='autoCrop':
- cyto=np.empty((len(CytoMaskChan),int(rect[2]),int(rect[3])),dtype=np.int16)
- for iChan in CytoMaskChan:
- cytoFull= tifffile.imread(imagePath, key=iChan)
- cyto[count,:,:] = cytoFull[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- count+=1
- else:
- cyto=np.empty((len(CytoMaskChan),rect[3],rect[2]),dtype=np.int16)
- for iChan in CytoMaskChan:
- cytoFull= tifffile.imread(imagePath, key=iChan)
- cyto[count,:,:] = cytoFull[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- count+=1
- cyto = np.amax(cyto,axis=0)
- cytoplasmMask,nucleiMaskTemp,cellMask = S3CytoplasmSegmentation(nucleiMask,cyto,TMAmask,args.cytoMethod,args.cytoDilation)
- exportMasks(nucleiMaskTemp,nucleiCrop,outputPath,filePrefix,'nuclei',commit,metadata,args.saveFig,args.saveMask)
- exportMasks(cytoplasmMask,cyto,outputPath,filePrefix,'cyto',commit,metadata,args.saveFig,args.saveMask)
- exportMasks(cellMask,cyto,outputPath,filePrefix,'cell',commit,metadata,args.saveFig,args.saveMask)
-
- cytoplasmMask,nucleiMaskTemp,cellMask = S3CytoplasmSegmentation(nucleiMask,cyto,TMAmask,'ring',args.cytoDilation)
- exportMasks(nucleiMaskTemp,nucleiCrop,outputPath,filePrefix,'nucleiRing',commit,metadata,args.saveFig,args.saveMask)
- exportMasks(cytoplasmMask,cyto,outputPath,filePrefix,'cytoRing',commit,metadata,args.saveFig,args.saveMask)
- exportMasks(cellMask,cyto,outputPath,filePrefix,'cellRing',commit,metadata,args.saveFig,args.saveMask)
-
- elif args.segmentCytoplasm == 'ignoreCytoplasm':
- exportMasks(nucleiMask,nucleiCrop,outputPath,filePrefix,'nuclei',commit,metadata)
- cellMask = nucleiMask
- exportMasks(nucleiMask,nucleiCrop,outputPath,filePrefix,'cell',commit,metadata)
- cytoplasmMask = nucleiMask
-
- detectPuncta = args.detectPuncta
- if (np.min(detectPuncta)>0):
- detectPuncta[:] = [number - 1 for number in detectPuncta] #convert 1-based indexing to 0-based indexing
- if len(detectPuncta) != len(args.punctaSigma):
- args.punctaSigma = args.punctaSigma[0] * np.ones(len(detectPuncta))
-
-
- if len(detectPuncta) != len(args.punctaSD):
- args.punctaSD = args.punctaSD[0] * np.ones(len(detectPuncta))
-
- counter=0
- for iPunctaChan in detectPuncta:
- punctaChan = tifffile.imread(imagePath,key = iPunctaChan)
- punctaChan = punctaChan[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3])]
- spots=S3punctaDetection(punctaChan,cellMask,args.punctaSigma[counter],args.punctaSD[counter])
- cellspotmask = nucleiMask
- P = regionprops(cellspotmask,intensity_image = spots ,cache=False)
- numSpots = []
- for prop in P:
- numSpots.append(np.uint16(np.round(prop.mean_intensity * prop.area)))
- np.savetxt(outputPath + os.path.sep + 'numSpots_chan' + str(iPunctaChan+1) +'.csv',(np.transpose(np.asarray(numSpots))),fmt='%10.5f')
- edges = 1-(cellMask>0)
- stacked_img=np.stack((np.uint16((spots+edges)>0)*np.amax(punctaChan),punctaChan),axis=0)
-
-
- outputPathPuncta = outputPath + os.path.sep + filePrefix + os.path.sep + 'punctaChan'+str(iPunctaChan+1) + 'Outlines.ome.tif'
-
- # metadata_args = dict(
- # pixel_sizes=(metadata.physical_size_y, metadata.physical_size_x),
- # pixel_size_units=('µm', 'µm'),
- # software= 's3segmenter v' + commit
- # )
- save_pyramid(
- stacked_img,
- outputPathPuncta,
- channel_names=['puncta outlines', 'image channel'],
- is_mask=False,
- **metadata
- )
-
- counter=counter+1
-
diff -r 41e8efe8df43 -r 96d0d969ebc9 macros.xml
--- a/macros.xml Fri Mar 11 23:37:49 2022 +0000
+++ b/macros.xml Fri Sep 16 20:05:54 2022 +0000
@@ -2,25 +2,36 @@
- labsyspharm/s3segmenter:@VERSION@
+
+ -->
+ labsyspharm/s3segmenter:@TOOL_VERSION@
-
- echo @VERSION@
-
- 1.3.12
- python3 $__tool_directory__/S3segmenter.py
+ 1.3.12
+ 0
+ 19.01
+
+
+ @CMD_BEGIN@ $S3SEG_CMD --help
+
diff -r 41e8efe8df43 -r 96d0d969ebc9 rowit.py
--- a/rowit.py Fri Mar 11 23:37:49 2022 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,76 +0,0 @@
-import numpy as np
-from skimage.util import view_as_windows, montage
-
-
-class WindowView(object):
-
- def __init__(
- self, img_shape, block_size, overlap_size
- ):
- self.img_shape = img_shape
- self.block_size = block_size
- self.overlap_size = overlap_size
-
- self.step_size = block_size - overlap_size
-
- def window_view_list(self, img, pad_mode='constant'):
- half = int(self.overlap_size / 2)
- img = np.pad(img, (
- (half, self.padded_shape[0] - self.img_shape[0] - half),
- (half, self.padded_shape[1] - self.img_shape[1] - half),
- ), mode=pad_mode)
-
- return self._window_view_list(img)
-
- def padding_mask(self):
- half = int(self.overlap_size / 2)
- padding_mask = np.ones(self.img_shape, dtype=np.bool)
- padding_mask = np.pad(padding_mask, (
- (half, self.padded_shape[0] - self.img_shape[0] - half),
- (half, self.padded_shape[1] - self.img_shape[1] - half),
- ), mode='constant', constant_values=0)
- return self._window_view_list(padding_mask)
-
- def reconstruct(self, img_window_view_list):
- grid_shape = self.window_view_shape[:2]
-
- start = int(self.overlap_size / 2)
- end = int(self.block_size - start)
-
- img_window_view_list = img_window_view_list[..., start:end, start:end]
-
- return montage(
- img_window_view_list, grid_shape=grid_shape
- )[:self.img_shape[0], :self.img_shape[1]]
-
- @property
- def padded_shape(self):
- padded_shape = np.array(self.img_shape) + self.overlap_size
- n = np.ceil((padded_shape - self.block_size) / self.step_size)
- padded_shape = (self.block_size + (n * self.step_size)).astype(np.int)
- return tuple(padded_shape)
-
- @property
- def window_view_shape(self):
- return view_as_windows(
- np.empty(self.padded_shape),
- self.block_size, self.step_size
- ).shape
-
- def _window_view_list(self, img):
- return (
- view_as_windows(img, self.block_size, self.step_size)
- .reshape(-1, self.block_size, self.block_size)
- )
-
-def crop_with_padding_mask(img, padding_mask, return_mask=False):
- if np.all(padding_mask == 1):
- return (img, padding_mask) if return_mask else img
- (r_s, r_e), (c_s, c_e) = [
- (i.min(), i.max() + 1)
- for i in np.where(padding_mask == 1)
- ]
- padded = np.zeros_like(img)
- img = img[r_s:r_e, c_s:c_e]
- padded[r_s:r_e, c_s:c_e] = 1
- return (img, padded) if return_mask else img
\ No newline at end of file
diff -r 41e8efe8df43 -r 96d0d969ebc9 s3segmenter.xml
--- a/s3segmenter.xml Fri Mar 11 23:37:49 2022 +0000
+++ b/s3segmenter.xml Fri Sep 16 20:05:54 2022 +0000
@@ -1,49 +1,49 @@
-
- S3segmenter is a Python-based set of functions that generates single cell (nuclei and cytoplasm) label masks.
+
+ single cell (nuclei and cytoplasm) label masks.
macros.xml
- @VERSION_CMD@
+
@@ -94,6 +87,9 @@
+
+
+
@@ -121,17 +117,24 @@
-
-
+
+
+
+
+
+
+
+
+
-
-
+
+
@@ -161,7 +164,29 @@
saveFig is True
+
+
+
+
+
+
+
+
+
diff -r 41e8efe8df43 -r 96d0d969ebc9 save_tifffile_pyramid.py
--- a/save_tifffile_pyramid.py Fri Mar 11 23:37:49 2022 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,114 +0,0 @@
-import numpy as np
-import tifffile
-import skimage.transform
-
-PHYSICAL_SIZE_UNIT = ['Ym', 'Zm', 'Em', 'Pm', 'Tm', 'Gm', 'Mm', 'km', 'hm', 'dam', 'm', 'dm', 'cm', 'mm', 'µm', 'nm', 'pm', 'fm', 'am', 'zm', 'ym', 'Å', 'thou', 'li', 'in', 'ft', 'yd', 'mi', 'ua', 'ly', 'pc', 'pt', 'pixel', 'reference frame']
-
-def normalize_image_shape(img):
- assert img.ndim in (2, 3), (
- 'image must be 2D (Y, X) or 3D (C, Y, X)'
- )
- if img.ndim == 2:
- img = img.reshape(1, *img.shape)
- assert np.argmin(img.shape) == 0, (
- '3D image must be in (C, Y, X) order'
- )
- return img
-
-def save_pyramid(
- out_img, output_path,
- pixel_sizes=(1, 1),
- pixel_size_units=('µm', 'µm'),
- channel_names=None,
- software=None,
- is_mask=False
-):
- assert '.ome.tif' in str(output_path)
- assert len(pixel_sizes) == len(pixel_size_units) == 2
- assert out_img.ndim in (2, 3), (
- 'image must be either 2D (Y, X) or 3D (C, Y, X)'
- )
-
- img_shape_ori = out_img.shape
- out_img = normalize_image_shape(out_img)
- img_shape = out_img.shape
-
- size_x, size_y = np.array(pixel_sizes, dtype=float)
- unit_x, unit_y = pixel_size_units
-
- assert (unit_x in PHYSICAL_SIZE_UNIT) & (unit_y in PHYSICAL_SIZE_UNIT), (
- f'pixel_size_units must be a tuple of the followings: '
- f'{", ".join(PHYSICAL_SIZE_UNIT)}'
- )
-
- n_channels = img_shape[0]
- if channel_names == None:
- channel_names = [f'Channel {i}' for i in range(n_channels)]
- else:
- if type(channel_names) == str:
- channel_names = [channel_names]
- n_channel_names = len(channel_names)
- assert n_channel_names == n_channels, (
- f'number of channel_names ({n_channel_names}) must match '
- f'number of channels ({n_channels})'
- )
-
- if software == None:
- software = ''
-
- metadata = {
- 'Creator': software,
- 'Pixels': {
- 'PhysicalSizeX': size_x,
- 'PhysicalSizeXUnit': unit_x,
- 'PhysicalSizeY': size_y,
- 'PhysicalSizeYUnit': unit_y,
- },
- 'Channel': {'Name': channel_names},
-
- }
-
- max_size = np.max(img_shape)
- subifds = np.ceil(np.log2(max_size / 1024)).astype(int)
-
- # use optimal tile size for disk space
- tile_size = 16*np.ceil(
- np.array(img_shape[1:]) / (2**subifds) / 16
- ).astype(int)
- options = {
- 'tile': tuple(tile_size)
- }
-
- with tifffile.TiffWriter(output_path, bigtiff=True) as tiff_out:
- tiff_out.write(
- data=out_img,
- metadata=metadata,
- software=software,
- subifds=subifds,
- **options
- )
- for i in range(subifds):
- if i == 0:
- down_2x_img = downsize_img_channels(out_img, is_mask=is_mask)
- else:
- down_2x_img = downsize_img_channels(down_2x_img, is_mask=is_mask)
- tiff_out.write(
- data=down_2x_img,
- subfiletype=1,
- **options
- )
-
- out_img = out_img.reshape(img_shape_ori)
- return
-
-def downsize_channel(img, is_mask):
- if is_mask:
- return img[::2, ::2]
- else:
- return skimage.transform.downscale_local_mean(img, (2, 2)).astype(img.dtype)
-
-def downsize_img_channels(img, is_mask):
- return np.array([
- downsize_channel(c, is_mask=is_mask)
- for c in img
- ])
diff -r 41e8efe8df43 -r 96d0d969ebc9 test-data/stack_probabilities.tiff
Binary file test-data/stack_probabilities.tiff has changed
diff -r 41e8efe8df43 -r 96d0d969ebc9 test-data/test.ome.tiff
Binary file test-data/test.ome.tiff has changed