changeset 0:37acf42a824b draft

"planemo upload for repository https://github.com/ohsu-comp-bio/S3segmenter commit d89a61efd4c465a1e6bf5b99b0f78fb19be5bdea-dirty"

--- /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

--- /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 @@
+<?xml version="1.0"?>
+<macros>
+    <xml name="requirements">
+        <requirements>
+            <requirement type="package" version="3.6">python</requirement>
+            <requirement type="package">scikit-learn</requirement>
+            <requirement type="package">scikit-image</requirement>
+            <requirement type="package">matplotlib</requirement>
+            <requirement type="package">tifffile</requirement>
+            <requirement type="package">opencv</requirement>
+        </requirements>
+    </xml>
+
+    <xml name="version_cmd">
+        <version_command>echo @VERSION@</version_command>
+    </xml>
+    <xml name="citations">
+        <citations>
+        </citations>
+    </xml>
+
+    <token name="@VERSION@">1.2.1</token>
+    <token name="@CMD_BEGIN@">python3 $__tool_directory__/S3segmenter.py</token>
+</macros>

--- /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

--- /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 @@
+<tool id="s3segmenter" name="s3segmenter" version="@VERSION@.3" profile="17.09">
+    <description>S3segmenter is a Python-based set of functions that generates single cell (nuclei and cytoplasm) label masks.</description>
+    <macros>
+        <import>macros.xml</import>
+    </macros>
+ 
+    <expand macro="requirements"/>
+    @VERSION_CMD@
+
+    <command detect_errors="exit_code"><![CDATA[
+
+        ln -s '${imagePath}' ./image.tif;
+
+        #if $contoursClassProbPath
+        ln -s ${contoursClassProbPath} ./ContoursPM.tif;
+        #end if
+
+        #if $nucleiClassProbPath
+        ln -s ${nucleiClassProbPath} ./NucleiPM.tif;
+        #end if
+
+        #if $stackProbPath
+        ln -s ${stackProbPath} ./Probabilities.tif;
+        #end if
+
+
+        @CMD_BEGIN@ 
+        --imagePath ./image.tif
+
+        #if $contoursClassProbPath
+        --contoursClassProbPath ./ContoursPM.tif
+        #end if
+
+        #if $nucleiClassProbPath
+        --nucleiClassProbPath ./NucleiPM.tif
+        #end if
+
+        #if $stackProbPath
+        --stackProbPath ./Probabilities.tif
+        #end if
+
+        --mask $mask
+        --probMapChan $probMapChan
+        --crop $crop
+        --cytoMethod $cytoMethod
+        --nucleiFilter $nucleiFilter
+        --nucleiRegion $nucleiRegion
+        --segmentCytoplasm $segmentCytoplasm
+        --cytoDilation $adv.cytoDilation
+        --logSigma $adv.logSigma
+        --CytoMaskChan $adv.CytoMaskChan
+        ##--TissueMaskChan $adv.TissueMaskChan
+        --detectPuncta $adv.detectPuncta
+        --punctaSigma $adv.punctaSigma
+        --punctaSD $adv.punctaSD
+
+        #if not $saveMask
+        --saveMask
+        #end if
+
+        #if not $saveFig
+        --saveFig
+        #end if
+
+        --outputPath .
+    ]]></command>
+
+
+    <inputs>
+
+        <param name="imagePath" type="data" format="tiff" label="Image File"/>
+        <param name="contoursClassProbPath" type="data" format="tiff" optional="true" label="Contours Class Probabilities"/>
+        <param name="nucleiClassProbPath" type="data" format="tiff" optional="true" label="Nuclei Class Probabilities"/>
+        <param name="stackProbPath" type="data" format="tiff" optional="true" label="Stack Probabilities"/>
+        <param name="mask" type="select" label="Choose mask: TMA, tissue, none.">
+            <option selected="true" value="tissue">tissue</option>
+            <option value="TMA">TMA</option>
+            <option value="none">none</option>
+        </param>
+        <param name="probMapChan" type="integer" value="-1" label="Probability Maps Channel"/>
+        <param name="crop" type="select" label="Crop Strategy">
+            <option selected="true" value="noCrop">No Crop</option>
+            <option value="autoCrop">Automatic Crop</option>
+            <option value="dearray">De-array</option>
+            <option value="plate">Plate</option>
+        </param>
+        <param name="cytoMethod" type="select" label="Cyto Method">
+            <option value="hybrid">Hybrid</option>
+            <option selected="true" value="distanceTransform">Distance Transform</option>
+            <option value="bwdistanceTransform">BW Distance Transform</option>
+            <option value="ring">Ring</option>
+        </param>
+        <param name="nucleiFilter" type="select" label="Nuclei Filter">
+            <option selected="true" value="IntPM">IntPM</option>
+            <option value="LoG">LoG</option>
+            <option value="Int">Int</option>
+            <option value="none">none</option>
+        </param>
+        <param name="nucleiRegion" type="select" label="Nuclei Region">
+            <option value="watershedContourDist">watershedContourDist</option>
+            <option selected="true" value="watershedContourInt">watershedContourInt</option>
+            <option value="watershedBWDist">watershedBWDist</option>
+            <option value="dilation">dilation</option>
+            <option value="localThreshold">localThreshold</option>
+        </param>
+        <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"/>
+
+        <section name="adv" title="Advanced Options" expanded="false">
+            <param name="cytoDilation" type="integer" value="5" label="Cyto Dilation"/>
+            <param name="logSigma" type="text" value="3 60" label="logSigma"/>
+            <param name="CytoMaskChan" type="text" value="1" label="Cyto Mask Channel"/>
+            <!-- Bug with S3Segmenter code, expects int not list
+            <param name="TissueMaskChan" type="text" value="-1" label="Tissue Mask Channel"/>
+            -->
+            <param name="detectPuncta" type="text" value="-1" label="Detect Puncta"/>
+            <param name="punctaSigma" type="text" value="1" label="Puncta Sigma"/>
+            <param name="punctaSD" type="text" value="4" label="Puncta SD"/>
+        </section>
+        
+
+    </inputs>
+    <outputs>
+        <data format="tiff" name="cell_mask" from_work_dir="*/cellMask.tif" label="cellMasks">
+            <filter>saveMask is True</filter>
+        </data>
+        <data format="tiff" name="nuclei_mask" from_work_dir="*/nucleiMask.tif" label="nucleiMasks">
+            <filter>saveMask is True</filter>
+        </data>
+        <data format="tiff" name="cell_outlines" from_work_dir="*/cellOutlines.tif" label="cellOutlines">
+            <filter>saveFig is True</filter>
+        </data>
+        <data format="tiff" name="nuclei_outlines" from_work_dir="*/nucleiOutlines.tif" label="nucleiOutlines">
+            <filter>saveFig is True</filter>
+        </data>
+    </outputs>
+    <help><![CDATA[
+Inputs are:
+
+1. an .ome.tif (preferably flat field corrected)
+2. a 3-class probability maps derived from a deep learning model such as UNet. Classes include background, nuclei contours, and nuclei foreground.
+
+The centers of each nuclei are obtained by finding local maxima from the nuclei foreground. These are used for marker-controlled watershed constrained by the nuclei contours.
+
+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>