diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b1d09ad..8f87339 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -48,10 +48,10 @@ jobs: # If the action is successful, the output will be available as a downloadable artifact - name: Upload processed result - uses: actions/upload-artifact@v2 + uses: actions/upload-artifact@v4 with: name: ex001-s3seg path: | ~/data/exemplar-001-cycle6/** ~/data/large/** - \ No newline at end of file + diff --git a/Dockerfile b/Dockerfile index ee7a03e..76df76a 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,6 +1,10 @@ -FROM python:3.8.0 +FROM python:3.13.9 -RUN pip install numpy==1.22.3 scikit-learn scikit-image==0.16.2 matplotlib tifffile==2021.6.6 opencv-python==4.3.0.36 ome_types==0.4.5 imagecodecs +RUN apt-get update \ + && DEBIAN_FRONTEND=noninteractive apt-get install -y libgl1 \ + && rm -rf /var/lib/apt/lists/* + +RUN pip install --no-cache numpy scikit-learn scikit-image matplotlib tifffile opencv-python ome_types imagecodecs COPY S3segmenter.py ./app/S3segmenter.py COPY save_tifffile_pyramid.py ./app/save_tifffile_pyramid.py diff --git a/S3segmenter.py b/S3segmenter.py index 351eb3a..3d5acac 100644 --- a/S3segmenter.py +++ b/S3segmenter.py @@ -8,6 +8,7 @@ from skimage.morphology import * from skimage.morphology import extrema from skimage import morphology +import skimage.measure._regionprops from skimage.measure import regionprops from skimage.transform import resize from skimage.filters import gaussian, threshold_otsu, threshold_local @@ -15,7 +16,7 @@ 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 scipy.ndimage import uniform_filter from os.path import * from os import listdir, makedirs, remove import pickle @@ -32,6 +33,7 @@ from save_tifffile_pyramid import save_pyramid import subprocess import ome_types +import threadpoolctl def imshowpair(A,B): @@ -58,6 +60,12 @@ def normI(I): J = J/(p99-p1); return J +def peak_local_max_mask(image, **kwargs): + idx = peak_local_max(image, **kwargs) + mask = np.zeros_like(image, dtype=bool) + mask[tuple(idx.T)] = True + return mask + def contour_pm_watershed( contour_pm, sigma=2, h=0, tissue_mask=None, padding_mask=None, min_area=None, max_area=None @@ -73,12 +81,11 @@ def contour_pm_watershed( tissue_mask, padding_mask ) - maxima = peak_local_max( + maxima = peak_local_max_mask( 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) @@ -99,10 +106,10 @@ def contour_pm_watershed( np.greater(maxima, 0, out=maxima) if padded is None: - return maxima.astype(np.bool) + return maxima.astype(bool) else: padded[padded == 1] = maxima.flatten() - return padded.astype(np.bool) + return padded.astype(bool) def S3AreaSegmenter(nucleiPM, images, TMAmask, threshold,fileprefix,outputPath): nucleiCenters = nucleiPM[:,:,0] @@ -152,7 +159,7 @@ def S3NucleiBypass(nucleiPM,nucleiImage,logSigma,TMAmask,nucleiFilter,nucleiRegi def props_of_keys(prop, keys): return [prop[k] for k in keys] - mean_ints, labels, areas = np.array(Parallel(n_jobs=1)(delayed(props_of_keys)(prop, prop_keys) + mean_ints, labels, areas = np.array(Parallel(n_jobs=1,verbose=joblibVerbose)(delayed(props_of_keys)(prop, prop_keys) for prop in P ) ).T @@ -171,7 +178,7 @@ def S3NucleiSegmentationWatershed(nucleiPM,nucleiImage,logSigma,TMAmask,nucleiFi 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 = peak_local_max_mask(extrema.h_maxima(255-nucleiContours,logSigma[0])) Imax = label(Imax).astype(np.int32) foregroundMask = watershed(nucleiContours, Imax, watershed_line=True) P = regionprops(foregroundMask, np.amax(nucleiCenters) - nucleiCenters - nucleiContours) @@ -179,7 +186,7 @@ def S3NucleiSegmentationWatershed(nucleiPM,nucleiImage,logSigma,TMAmask,nucleiFi 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) + mean_ints, labels, areas = np.array(Parallel(n_jobs=numWorkers,verbose=joblibVerbose)(delayed(props_of_keys)(prop, prop_keys) for prop in P ) ).T @@ -210,9 +217,10 @@ def props_of_keys(prop, keys): maxArea = (logSigma[1]**2)*3/4 minArea = (logSigma[0]**2)*3/4 - + + print(datetime.datetime.now(), '==> Performing watershed') foregroundMask = np.array( - Parallel(n_jobs=6)(delayed(contour_pm_watershed)( + Parallel(n_jobs=numWorkers,verbose=joblibVerbose)(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)) @@ -228,7 +236,7 @@ def props_of_keys(prop, keys): elif nucleiFilter == 'Int': int_img = nucleiImage - print(' ', datetime.datetime.now(), 'regionprops') + print(datetime.datetime.now(), '==> Filtering out low-confidence events') P = regionprops(foregroundMask, int_img) def props_of_keys(prop, keys): @@ -236,7 +244,7 @@ def props_of_keys(prop, 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) + Parallel(n_jobs=numWorkers,verbose=joblibVerbose)(delayed(props_of_keys)(prop, prop_keys) for prop in P ) ).T @@ -349,7 +357,7 @@ def exportMasks(mask,image,outputPath,filePrefix,fileName,commit,metadata_args,s 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 + fgm=peak_local_max_mask(Ilog, footprint=np.ones((3, 3)))*tissueMask test=Ilog[fgm==1] med = np.median(test) mad = np.median(np.absolute(test - med)) @@ -385,15 +393,32 @@ def S3punctaDetection(spotChan,mask,sigma,SD): parser.add_argument("--punctaSD", nargs = '+', type=float, default=[4]) parser.add_argument("--saveMask",action='store_false') parser.add_argument("--saveFig",action='store_false') + parser.add_argument("--numWorkers", type=int, default=0, help="Parallel worker processes, recommend <= 4 (default=autodetect)") + parser.add_argument("--verbose", action="store_true", help="Enable for extra logging output") args = parser.parse_args() - + imagePath = args.imagePath outputPath = args.outputPath nucleiClassProbPath = args.nucleiClassProbPath contoursClassProbPath = args.contoursClassProbPath stackProbPath = args.stackProbPath maskPath = args.maskPath - + joblibVerbose = 10 if args.verbose else 0 + + threadpoolctl.threadpool_limits(1) + if args.numWorkers == 0: + if hasattr(os, 'sched_getaffinity'): + numCpus = len(os.sched_getaffinity(0)) + else: + numCpus = multiprocessing.cpu_count() + numWorkers = min(numCpus, 4) + print( + f"Using {numWorkers} worker processes based on available CPU count ({numCpus})." + ) + print() + else: + numWorkers = args.numWorkers + commit = '1.3.11'#subprocess.check_output(['git', 'describe', '--tags']).decode('ascii').strip() metadata = getMetadata(imagePath,commit) @@ -438,6 +463,8 @@ def S3punctaDetection(spotChan,mask,sigma,SD): TissueMaskChan[:] = [number - 1 for number in TissueMaskChan]#convert 1-based indexing to 0-based indexing TissueMaskChan.append(nucMaskChan) + print(datetime.datetime.now(), '==> Preprocessing images') + #crop images if needed if args.crop == 'interactiveCrop': nucleiCrop = tifffile.imread(imagePath,key = nucMaskChan) @@ -513,6 +540,7 @@ def S3punctaDetection(spotChan,mask,sigma,SD): del nucleiPM # cytoplasm segmentation if args.segmentCytoplasm == 'segmentCytoplasm': + print(datetime.datetime.now(), '==> Segmenting cytoplasm') 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) @@ -550,6 +578,7 @@ def S3punctaDetection(spotChan,mask,sigma,SD): detectPuncta = args.detectPuncta if (np.min(detectPuncta)>0): + print(datetime.datetime.now(), '==> Detecting puncta') 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)) @@ -606,4 +635,6 @@ def S3punctaDetection(spotChan,mask,sigma,SD): channel_names= spot_channel_names, is_mask=False, **metadata - ) \ No newline at end of file + ) + + print(datetime.datetime.now(), '==> Done!') diff --git a/__pycache__/rowit.cpython-36.pyc b/__pycache__/rowit.cpython-36.pyc deleted file mode 100644 index 5a2c4f3..0000000 Binary files a/__pycache__/rowit.cpython-36.pyc and /dev/null differ diff --git a/__pycache__/rowit.cpython-38.pyc b/__pycache__/rowit.cpython-38.pyc deleted file mode 100644 index 76a124f..0000000 Binary files a/__pycache__/rowit.cpython-38.pyc and /dev/null differ diff --git a/rowit.py b/rowit.py index 2e9ab1c..14e0482 100644 --- a/rowit.py +++ b/rowit.py @@ -47,7 +47,7 @@ def reconstruct(self, img_window_view_list): 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) + padded_shape = (self.block_size + (n * self.step_size)).astype(int) return tuple(padded_shape) @property @@ -73,4 +73,4 @@ def crop_with_padding_mask(img, padding_mask, return_mask=False): 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 + return (img, padded) if return_mask else img