Mercurial > repos > perssond > s3segmenter
changeset 2:96d0d969ebc9 draft default tip
planemo upload for repository https://github.com/goeckslab/tools-mti/tree/main/tools/s3segmenter commit 0f4f17235c5961c2fd3d4c30180507f66214c11d
| author | goeckslab |
|---|---|
| date | Fri, 16 Sep 2022 20:05:54 +0000 |
| parents | 41e8efe8df43 |
| children | |
| files | S3segmenter.py macros.xml rowit.py s3segmenter.xml save_tifffile_pyramid.py test-data/stack_probabilities.tiff test-data/test.ome.tiff |
| diffstat | 7 files changed, 70 insertions(+), 817 deletions(-) [+] |
line wrap: on
line diff
--- 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 -
--- 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 @@ <macros> <xml name="requirements"> <requirements> - <container type="docker">labsyspharm/s3segmenter:@VERSION@</container> + <!-- <requirement type="package" version="3.7">python</requirement> <requirement type="package">scikit-learn</requirement> <requirement version="0.14.2" type="package">scikit-image</requirement> <requirement type="package">matplotlib</requirement> <requirement version="2021.6.6" type="package">tifffile</requirement> <requirement type="package">opencv</requirement> - <!-- <requirement type="package">ome_types</requirement> --> + --> + <container type="docker">labsyspharm/s3segmenter:@TOOL_VERSION@</container> </requirements> </xml> - <xml name="version_cmd"> - <version_command>echo @VERSION@</version_command> - </xml> <xml name="citations"> <citations> </citations> </xml> - <token name="@VERSION@">1.3.12</token> - <token name="@CMD_BEGIN@">python3 $__tool_directory__/S3segmenter.py</token> + <token name="@TOOL_VERSION@">1.3.12</token> + <token name="@VERSION_SUFFIX@">0</token> + <token name="@PROFILE@">19.01</token> + <token name="@CMD_BEGIN@"><![CDATA[ + S3SEG_CMD="" && + if [ -f "/app/S3segmenter.py" ]; then + export S3SEG_CMD="python /app/S3segmenter.py"; + export PYTHONPATH="/app"; + else + export S3SEG_CMD="S3segmenter.py"; + fi && + ]]></token> + <xml name="version_cmd"> + <version_command>@CMD_BEGIN@ $S3SEG_CMD --help</version_command> + </xml> </macros>
--- 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
--- 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 @@ -<tool id="s3segmenter" name="s3segmenter" version="@VERSION@.0" profile="17.09"> - <description>S3segmenter is a Python-based set of functions that generates single cell (nuclei and cytoplasm) label masks.</description> +<tool id="s3segmenter" name="s3segmenter" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="@PROFILE@"> + <description>single cell (nuclei and cytoplasm) label masks.</description> <macros> <import>macros.xml</import> </macros> <expand macro="requirements"/> - @VERSION_CMD@ + <expand macro="version_cmd"/> <command detect_errors="exit_code"><![CDATA[ - ln -s '${imagePath}' ./image.tif; + ln -s '${imagePath}' './image.tif' && #if $contoursClassProbPath - ln -s ${contoursClassProbPath} ./ContoursPM.tif; + ln -s '${contoursClassProbPath}' './ContoursPM.tif' && #end if #if $nucleiClassProbPath - ln -s ${nucleiClassProbPath} ./NucleiPM.tif; + ln -s '${nucleiClassProbPath}' './NucleiPM.tif' && #end if #if $stackProbPath - ln -s ${stackProbPath} ./Probabilities.tif; + ln -s '${stackProbPath}' './Probabilities.tif' && #end if - @CMD_BEGIN@ - --imagePath ./image.tif + @CMD_BEGIN@ \$S3SEG_CMD + --imagePath './image.tif' #if $contoursClassProbPath - --contoursClassProbPath ./ContoursPM.tif + --contoursClassProbPath './ContoursPM.tif' #end if #if $nucleiClassProbPath - --nucleiClassProbPath ./NucleiPM.tif + --nucleiClassProbPath './NucleiPM.tif' #end if #if $stackProbPath - --stackProbPath ./Probabilities.tif + --stackProbPath './Probabilities.tif' #end if --probMapChan $probMapChan --crop $crop_select.crop #if $crop_select.crop == "dearray" - --maskPath $crop_select.maskPath + --maskPath '$crop_select.maskPath' #end if --cytoMethod $cytoMethod @@ -63,16 +63,9 @@ --detectPuncta $adv.detectPuncta --punctaSigma $adv.punctaSigma --punctaSD $adv.punctaSD - - #if not $saveMask - --saveMask - #end if - - #if not $saveFig - --saveFig - #end if - - --outputPath . + $saveMask + $saveFig + --outputPath '.' ]]></command> @@ -94,6 +87,9 @@ <when value="dearray"> <param name="maskPath" type="data" format="tiff" label="TMA Mask File"/> </when> + <when value="noCrop" /> + <when value="autoCrop" /> + <when value="plate" /> </conditional> <param name="cytoMethod" type="select" label="Cyto Method"> @@ -121,17 +117,24 @@ <option value="pixellevel">pixellevel</option> </param> <when value="pixellevel"> - <param name="pixelThreshold" type="float" value="-1.0" Label="Pixel Threshold"/> - <param name="pixelMaskChan" type="text" value="2" Label="Pixel Mask Channel"/> + <param name="pixelThreshold" type="float" value="-1.0" label="Pixel Threshold"/> + <param name="pixelMaskChan" type="text" value="2" label="Pixel Mask Channel"/> </when> + <when value="watershedContourDist"/> + <when value="watershedContourInt"/> + <when value="watershedBWDist"/> + <when value="dilation"/> + <when value="localThreshold"/> + <when value="localMax"/> + <when value="bypass"/> </conditional> <param name="segmentCytoplasm" type="select" label="Segment Cytoplasm"> <option value="segmentCytoplasm">segmentCytoplasm</option> <option selected="true" value="ignoreCytoplasm">ignoreCytoplasm</option> </param> - <param name="saveMask" type="boolean" checked="true" label="Save Mask"/> - <param name="saveFig" type="boolean" checked="true" label="Save Figure"/> + <param argument="saveMask" type="boolean" checked="true" truevalue="--saveMask" falsevalue="" label="Save Mask"/> + <param argument="saveFig" type="boolean" checked="true" truevalue="--saveFig" falsevalue="" label="Save Figure"/> <section name="adv" title="Advanced Options" expanded="false"> <param name="cytoDilation" type="integer" value="5" label="Cyto Dilation"/> @@ -161,7 +164,29 @@ <filter>saveFig is True</filter> </data> </outputs> + <tests> + <test> + <param name="imagePath" value="test.ome.tiff" /> + <param name="stackProbPath" value="stack_probabilities.tiff" /> + <param name="punctaSD" value="4" /> + <output name="cell_mask" ftype="tiff"> + <assert_contents> + <has_size value="6600000" delta="1000000" /> + </assert_contents> + </output> + <output name="nuclei_mask" ftype="tiff"> + <assert_contents> + <has_size value="6600000" delta="1000000" /> + </assert_contents> + </output> + </test> + </tests> <help><![CDATA[ +------------------- +S3segmenter +------------------- +**S3segmenter** is a Python-based set of functions that generates single cell (nuclei and cytoplasm) label masks. + Inputs are: 1. an .ome.tif (preferably flat field corrected) @@ -172,7 +197,6 @@ To segment cytoplasm, the nuclei are in turn used for a marker-controlled watershed segmentation constrained by a cytoplasmic marker such as B-catenin. The channel number of this marker must be specified. A 3-pixel annulus around each nucleus will also be used to segment cytoplasm. The source repository can be found here: https://github.com/HMS-IDAC/S3segmenter -OHSU Wrapper Repo: https://github.com/ohsu-comp-bio/S3segmenter ]]></help> <expand macro="citations" /> </tool>
--- 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 - ])