# HG changeset patch # User perssond # Date 1615508320 0 # Node ID 37acf42a824b20cd0203330f9cde8629bedae317 "planemo upload for repository https://github.com/ohsu-comp-bio/S3segmenter commit d89a61efd4c465a1e6bf5b99b0f78fb19be5bdea-dirty" diff -r 000000000000 -r 37acf42a824b S3segmenter.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/S3segmenter.py Fri Mar 12 00:18:40 2021 +0000 @@ -0,0 +1,629 @@ +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 + + +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,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))) + np.savetxt(outputPath + os.path.sep + 'area.csv',(np.transpose(np.asarray(area))),fmt='%10.5f') + return TMAmask + + +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': + 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) + elif nucleiRegion =='bypass': + foregroundMask = nucleiPM[:,:,0] + 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) + 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,saveFig=True,saveMasks = True): + outputPath =outputPath + os.path.sep + filePrefix + if not os.path.exists(outputPath): + os.makedirs(outputPath) + if saveMasks ==True: + kwargs={} + kwargs['bigtiff'] = True + kwargs['photometric'] = 'minisblack' + resolution = np.round(1) + kwargs['resolution'] = (resolution, resolution, 'cm') + kwargs['metadata'] = None + kwargs['description'] = '!!xml!!' + imsave(outputPath + os.path.sep + fileName + 'Mask.tif',mask, plugin="tifffile") + + 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) + tifffile.imsave(outputPath + os.path.sep + fileName + 'Outlines.tif',stacked_img) + +def S3punctaDetection(spotChan,mask,sigma,SD): + Ilog = -gaussian_laplace(np.float32(spotChan),sigma) + fgm=peak_local_max(Ilog, indices=False,footprint=np.ones((3, 3))) + 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) + +# stacked_img=np.stack((spots*4095,nucleiCrop),axis=0) +# tifffile.imsave('D:/Seidman/ZeissTest Sets/registration/spottest.tif',stacked_img) + + + # assign nan to tissue mask + + +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','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=[1]) + parser.add_argument("--pixelMaskChan",type=int, nargs = '+', default=[1]) + parser.add_argument("--TissueMaskChan",type=int, nargs = '+', default=-1) + parser.add_argument("--detectPuncta",type=int, nargs = '+', default=[-1]) + parser.add_argument("--punctaSigma", nargs = '+', type=float, default=[1]) + 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() + + # gather filename information + #exemplar001 +# imagePath = 'D:/LSP/cycif/testsets/exemplar-001/registration/exemplar-001small.ome.tif' +# outputPath = 'D:/LSP/cycif/testsets/exemplar-001/segmentation' +## nucleiClassProbPath = 'D:/LSP/cycif/testsets/exemplar-001/probability_maps/exemplar-001_NucleiPM_0.tif' +## contoursClassProbPath = 'D:/LSP/cycif/testsets/exemplar-001/probability_maps/exemplar-001_ContoursPM_0.tif' +# contoursClassProbPath ='' +# stackProbPath = 'D:/LSP/cycif/testsets/exemplar-001_Probabilities_0.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-001/dearray/masks/A1_mask.tif' +# args.cytoMethod = 'hybrid' +# args.nucleiRegion = 'localThreshold' +# args.segmentCytoplasm = 'segmentCytoplasm' + +# exemplar002 +# imagePath = 'D:/LSP/cycif/testsets/exemplar-002/dearray/1.tif' +# outputPath = 'D:/LSP/cycif/testsets/exemplar-002/segmentation' +# nucleiClassProbPath = ''#'D:/LSP/cycif/testsets/exemplar-002/prob_map/1_NucleiPM_1.tif' +# contoursClassProbPath = ''#'D:/LSP/cycif/testsets/exemplar-002/prob_map/1_ContoursPM_1.tif' +# stackProbPath = 'D:/LSP/cycif/testsets/exemplar-002/probability-maps/unmicst_1new_Probabilities_1.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-002/dearrayPython/masks/1_mask.tif' +# args.probMapChan =1 +# args.cytoMethod = 'hybrid' +# args.mask = 'TMA' +# args.crop = 'dearray' +# args.crop = 'autoCrop' +# args.segmentCytoplasm = 'segmentCytoplasm' + #PTCL +# imagePath = 'D:/LSP/cycif/testsets/PTCL/dearrayPython/1.tif' +# outputPath = 'D:/LSP/cycif/testsets/PTCL/dearrayPython/segmentation' +# nucleiClassProbPath = ''#'D:/LSP/cycif/testsets/exemplar-002/prob_map/1_NucleiPM_1.tif' +# contoursClassProbPath = ''#'D:/LSP/cycif/testsets/exemplar-002/prob_map/1_ContoursPM_1.tif' +# stackProbPath = 'D:/LSP/cycif/testsets/PTCL/prob_maps/1_Probabilities_40.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-002/dearrayPython/masks/1_mask.tif' +# args.crop = 'dearray' +# args.segmentCytoplasm = 'segmentCytoplasm' +# + #punctatest +# imagePath = 'Z:/IDAC/Clarence/Seidman/ZeissMouseTestSet/31082020_XXWT_Txnip550_Ddx3x647_WGA488_40x_1.tif' +# outputPath = 'Z:/IDAC/Clarence/Seidman/ZeissMouseTestSet/segmentation' +## nucleiClassProbPath = 'Z:/IDAC/Clarence/Seidman/ZeissMouseTestSet/probability-maps/test_NucleiPM_1.tif' +## contoursClassProbPath = 'Z:/IDAC/Clarence/Seidman/DanMouse/probability-maps/test_ContoursPM_1.tif' +# contoursClassProbPath ='' +# stackProbPath = 'Z:/IDAC/Clarence/Seidman/ZeissMouseTestSet/probability-maps/31082020_XXWT_Txnip550_Ddx3x647_WGA488_40x_1_Probabilities_2.tif' +# maskPath = 'D:/Seidman/ZeissTest Sets/segmentation/13042020_15AP_FAP488_LINC550_DCN647_WGA_40x_1/cellMask.tif' +# args.nucleiRegion = 'localThreshold' +# args.logSigma = [45, 300] +# args.segmentCytoplasm = 'ignoreCytoplasm' +# args.detectPuncta = [1] +# args.punctaSigma = [1] +# args.punctaSD = [3] + + + #plate +# imagePath = 'Y:/sorger/data/computation/Jeremy/caitlin-ddd-cycif-registered/Plate1/E3_fld_1/registration/E3_fld_1.ome.tif' +# outputPath = 'Y:/sorger/data/computation/Jeremy/caitlin-ddd-cycif-registered/Plate1/E3_fld_1/segmentation' +# nucleiClassProbPath = 'Y:/sorger/data/computation/Jeremy/caitlin-ddd-cycif-registered/Plate1/E3_fld_1/prob_maps/E3_fld_1_NucleiPM_1.tif' +# contoursClassProbPath = 'Y:/sorger/data/computation/Jeremy/caitlin-ddd-cycif-registered/Plate1/E3_fld_1/prob_maps/E3_fld_1_ContoursPM_1.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-001/dearray/masks/A1_mask.tif' +# args.crop = 'plate' +# args.cytoMethod ='hybrid' + + #large tissue +# imagePath = 'D:/WD-76845-097.ome.tif' +# outputPath = 'D:/' +## nucleiClassProbPath = 'D:/WD-76845-097_NucleiPM_28.tif' +# contoursClassProbPath = ""#'D:/WD-76845-097_ContoursPM_28.tif' +# stackProbPath = 'D:/ilastik/WD-76845-097_Probabilities_0.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-001/dearray/masks/A1_mask.tif' +# args.nucleiRegion = 'localThreshold' +# args.crop = 'autoCrop' +# args.TissueMaskChan =[0] + + #maskRCNN +# imagePath = 'D:/Seidman/s3segtest/registration/1.tif' +# outputPath = 'D:/Seidman/s3segtest/segmentation' +# stackProbPath = 'D:/Seidman/s3segtest/1_Probabilities_0.tif' +# contoursClassProbPath = '' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-001/dearray/masks/A1_mask.tif' +# args.nucleiRegion = 'bypass' +# args.crop = 'noCrop' +# args.TissueMaskChan =[0] +# args.logSigma = [45,300] + +# #pixellevel +# imagePath = 'D:/Olesja/D2107_spleen_DAPI-EdU_01/Layer0/D2107_spleen_DAPI-EdU_01.btf' +# outputPath = 'D:/Olesja/D2107_spleen_DAPI-EdU_01/segmentation' +# stackProbPath = 'D:/Seidman/s3segtest/1_Probabilities_0.tif' +# contoursClassProbPath = 'D:/Olesja/D2107_spleen_DAPI-EdU_01/prob_maps3Class/D2107_spleen_DAPI-EdU_01_ContoursPM_0.tif' +# nucleiClassProbPath = 'D:/Olesja/D2107_spleen_DAPI-EdU_01/prob_maps3Class/D2107_spleen_DAPI-EdU_01_NucleiPM_0.tif' +# maskPath = 'D:/LSP/cycif/testsets/exemplar-001/dearray/masks/A1_mask.tif' +# args.nucleiRegion = 'pixellevel' +# args.crop = 'noCrop' +# args.TissueMaskChan =[0] +# args.pixelThreshold = 0.06 + + imagePath = args.imagePath + outputPath = args.outputPath + nucleiClassProbPath = args.nucleiClassProbPath + contoursClassProbPath = args.contoursClassProbPath + stackProbPath = args.stackProbPath + maskPath = args.maskPath + + fileName = os.path.basename(imagePath) + filePrefix = fileName[0:fileName.index('.')] + + 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) + index = re.search('%s(.*)%s' % ('Probabilities', '.tif'), stackProbPath).group(1) + if len(index)==0: + nucMaskChan = args.probMapChan + else: + nucMaskChan = int(re.sub("\D", "", index)) + else: + print('NO PROBABILITY MAP PROVIDED') + if args.probMapChan ==-1: + if nucMaskChan ==-1: + sys.exit('INVALID NUCLEI CHANNEL SELECTION. SELECT CHANNEL USING --probMapChan') + else: + print('extracting nuclei channel from filename') + else: + nucMaskChan = args.probMapChan + + if args.TissueMaskChan==-1: + TissueMaskChan = copy.copy(args.CytoMaskChan) + TissueMaskChan.append(nucMaskChan) + else: + TissueMaskChan = args.TissueMaskChan[:] + 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) + nucleiPM = nucleiProbMaps[int(PMrect[0]):int(PMrect[0]+PMrect[2]), int(PMrect[1]):int(PMrect[1]+PMrect[3]),0:2] + 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) + + # tissue_gauss = tissueCrop +# tissue_gauss1 = tissue_gauss.astype(float) +# tissue_gauss1[tissue_gauss>np.percentile(tissue_gauss,99)]=np.nan +# TMAmask = np.log2(tissue_gauss+1)>threshold_otsu(np.log2(tissue_gauss+1)) + #imshow(TMAmask) + del tissue_gauss, tissue + + # nuclei segmentation + if args.nucleiRegion == 'pixellevel': + pixelTissue = np.empty((len(args.pixelMaskChan),nucleiCrop.shape[0],nucleiCrop.shape[1]),dtype=np.uint16) + count=0 + for iChan in args.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,outputPath) + 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(args.CytoMaskChan),nucleiCrop.shape[0],nucleiCrop.shape[1]),dtype=np.uint16) + for iChan in args.CytoMaskChan: + cyto[count,:,:] = tifffile.imread(imagePath, key=iChan) + count+=1 + elif args.crop =='autoCrop': + cyto=np.empty((len(args.CytoMaskChan),int(rect[2]),int(rect[3])),dtype=np.int16) + for iChan in args.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(args.CytoMaskChan),rect[3],rect[2]),dtype=np.int16) + for iChan in args.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',args.saveFig,args.saveMask) + exportMasks(cytoplasmMask,cyto,outputPath,filePrefix,'cyto',args.saveFig,args.saveMask) + exportMasks(cellMask,cyto,outputPath,filePrefix,'cell',args.saveFig,args.saveMask) + + cytoplasmMask,nucleiMaskTemp,cellMask = S3CytoplasmSegmentation(nucleiMask,cyto,TMAmask,'ring',args.cytoDilation) + exportMasks(nucleiMaskTemp,nucleiCrop,outputPath,filePrefix,'nucleiRing',args.saveFig,args.saveMask) + exportMasks(cytoplasmMask,cyto,outputPath,filePrefix,'cytoRing',args.saveFig,args.saveMask) + exportMasks(cellMask,cyto,outputPath,filePrefix,'cellRing',args.saveFig,args.saveMask) + + elif args.segmentCytoplasm == 'ignoreCytoplasm': + exportMasks(nucleiMask,nucleiCrop,outputPath,filePrefix,'nuclei') + cellMask = nucleiMask + exportMasks(nucleiMask,nucleiCrop,outputPath,filePrefix,'cell') + cytoplasmMask = nucleiMask + + + if (np.min(args.detectPuncta)>-1): + if len(args.detectPuncta) != len(args.punctaSigma): + args.punctaSigma = args.punctaSigma[0] * np.ones(len(args.detectPuncta)) + + + if len(args.detectPuncta) != len(args.punctaSD): + args.punctaSD = args.punctaSD[0] * np.ones(len(args.detectPuncta)) + + counter=0 + for iPunctaChan in args.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#tifffile.imread(maskPath) #REMOVE THIS LATER + 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) +'.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) + tifffile.imsave(outputPath + os.path.sep + filePrefix + os.path.sep + 'punctaChan'+str(iPunctaChan) + 'Outlines.tif',stacked_img) + counter=counter+1 + #fix bwdistance watershed diff -r 000000000000 -r 37acf42a824b macros.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros.xml Fri Mar 12 00:18:40 2021 +0000 @@ -0,0 +1,24 @@ + + + + + python + scikit-learn + scikit-image + matplotlib + tifffile + opencv + + + + + echo @VERSION@ + + + + + + + 1.2.1 + python3 $__tool_directory__/S3segmenter.py + diff -r 000000000000 -r 37acf42a824b rowit.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rowit.py Fri Mar 12 00:18:40 2021 +0000 @@ -0,0 +1,76 @@ +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 000000000000 -r 37acf42a824b s3segmenter.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/s3segmenter.xml Fri Mar 12 00:18:40 2021 +0000 @@ -0,0 +1,155 @@ + + S3segmenter is a Python-based set of functions that generates single cell (nuclei and cytoplasm) label masks. + + macros.xml + + + + @VERSION_CMD@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + +
+ + + + + + saveMask is True + + + saveMask is True + + + saveFig is True + + + saveFig is True + + + + +