diff --git a/README.md b/README.md index 1396980..992bf7d 100644 --- a/README.md +++ b/README.md @@ -116,7 +116,7 @@ Running the following command to generate the configuration file palom-svs-helper -i "Y:\DATA\SARDANA\MIHC\768473\RAW" -n "*Hem*" -o "Y:\DATA\SARDANA\MIHC\768473\RAW\palom\768473.ome.tif" -c "Y:\DATA\SARDANA\MIHC\768473\768473.yml" ``` -And the resulting `Y:\DATA\SARDANA\MIHC\768473\768473.yml` file +And the resulting `Y:\DATA\SARDANA\MIHC\768473\768473.yml` file ```yaml input dir: Y:\DATA\SARDANA\MIHC\768473\RAW @@ -153,7 +153,7 @@ palom-svs run -c "Y:\DATA\SARDANA\MIHC\768473\768473.yml" ``` When the process is finished, a pyramidal OME-TIFF file will be generated along -with PNG files showing the feature-based registration results and a log file. +with PNG files showing the feature-based registration results and a log file. ``` Y:\DATA\SARDANA\MIHC\768473\RAW @@ -199,7 +199,7 @@ c1rp = palom.color.PyramidHaxProcessor(c1r.pyramid, thumbnail_level=THUMBNAIL_LE c2rp = palom.color.PyramidHaxProcessor(c2r.pyramid, thumbnail_level=THUMBNAIL_LEVEL) c21l = palom.align.Aligner( - c1rp.get_processed_color(LEVEL), + c1rp.get_processed_color(LEVEL), c2rp.get_processed_color(LEVEL), ref_thumbnail=c1rp.get_processed_color(THUMBNAIL_LEVEL).compute(), moving_thumbnail=c2rp.get_processed_color(THUMBNAIL_LEVEL).compute(), diff --git a/palom/__init__.py b/palom/__init__.py index 66dfb6d..52c2b59 100644 --- a/palom/__init__.py +++ b/palom/__init__.py @@ -11,7 +11,7 @@ align, pyramid, color, - + # debugging block_affine, extract, diff --git a/palom/align.py b/palom/align.py index 2750582..5a9685c 100644 --- a/palom/align.py +++ b/palom/align.py @@ -93,12 +93,12 @@ def viz_shifts(shifts, grid_shape, dcenter=None, ax=None): cax = ax.inset_axes([1.04, 0.0, 0.02, 1]) colorbar = plt.colorbar(im, cax=cax) colorbar.set_ticks(colorbar_ticks) - + return ax def block_affine_matrices(mx, shifts): - + def shift_affine_mx(mx, shift): y, x = shift mx_shift = np.eye(3) @@ -155,7 +155,7 @@ def coarse_register_affine(self, **kwargs): self.coarse_affine_matrix = np.vstack( [affine_matrix, [0, 0, 1]] ) - + @property def affine_matrix(self): if not hasattr(self, 'coarse_affine_matrix'): @@ -169,13 +169,13 @@ def affine_matrix(self): mx_moving ) return affine_matrix - + @property def tform(self): return skimage.transform.AffineTransform( matrix=self.affine_matrix ) - + def affine_transformed_moving_img(self, mxs=None): if mxs is None: mxs = self.affine_matrix @@ -185,7 +185,7 @@ def affine_transformed_moving_img(self, mxs=None): return block_affine_transformed_moving_img( ref_img, moving_img, mxs ) - + def compute_shifts(self, pcc_kwargs=None): logger.info(f"Computing block-wise shifts") ref_img = self.ref_img @@ -212,7 +212,7 @@ def constrain_shifts(self): self.original_shifts, self.grid_shape ) - + @property def block_affine_matrices(self): mx = self.affine_matrix @@ -226,7 +226,7 @@ def block_affine_matrices_da(self): self.grid_shape ) - def overlay_grid(self, ax=None): + def overlay_grid(self, ax=None): import matplotlib.pyplot as plt img = self.ref_thumbnail img = skimage.exposure.rescale_intensity(img, out_range=np.uint16) @@ -249,7 +249,7 @@ def overlay_grid(self, ax=None): # checkerboard pattern ax.imshow(np.indices(shape).sum(axis=0) % 2, cmap='cool', alpha=0.2) return grid - + def plot_shifts(self): import matplotlib.pyplot as plt fig, axs = plt.subplots(1, 2, sharex=True, sharey=True) @@ -271,14 +271,14 @@ def get_aligner( ) if thumbnail_level1 <= -1: thumbnail_level1 += len(reader1.pyramid) if thumbnail_level2 <= -1: thumbnail_level2 += len(reader2.pyramid) - return Aligner( - reader1.read_level_channels(level1, channel1), - reader2.read_level_channels(level2, channel2), - reader1.read_level_channels(thumbnail_level1, channel1), - reader2.read_level_channels(thumbnail_level2, channel2), - reader1.level_downsamples[thumbnail_level1] / reader1.level_downsamples[level1], - reader2.level_downsamples[thumbnail_level2] / reader2.level_downsamples[level2] - ) + return Aligner( + reader1.read_level_channels(level1, channel1), + reader2.read_level_channels(level2, channel2), + reader1.read_level_channels(thumbnail_level1, channel1), + reader2.read_level_channels(thumbnail_level2, channel2), + reader1.level_downsamples[thumbnail_level1] / reader1.level_downsamples[level1], + reader2.level_downsamples[thumbnail_level2] / reader2.level_downsamples[level2] + ) def match_thumbnail_level(readers): diff --git a/palom/block_affine.py b/palom/block_affine.py index d074065..5a999cb 100644 --- a/palom/block_affine.py +++ b/palom/block_affine.py @@ -30,13 +30,13 @@ def block_affine_dask( def block_affine( - position, block_shape, + position, block_shape, transformation, src_img, fill_empty=0, multichannel=False, is_mask=False ): assert np.min(block_shape) >= 0, ( - f'block_shape {block_shape} is invalid' + f'block_shape {block_shape} is invalid' ) if multichannel: assert np.min(block_shape) == block_shape[0], ( @@ -68,7 +68,7 @@ def block_affine( x0_src, y0_src = np.floor( np.clip(inversed_corners.min(axis=0), 0, None) ).astype(int) - + if multichannel: src_img_block = src_img[:, y0_src:y1_src, x0_src:x1_src] src_img_block = np.moveaxis(src_img_block, 0, 2) diff --git a/palom/cli/align_he.py b/palom/cli/align_he.py index 2262f69..76d1576 100644 --- a/palom/cli/align_he.py +++ b/palom/cli/align_he.py @@ -42,16 +42,16 @@ def align_he( aligner.coarse_register_affine(n_keypoints=5000, detect_flip_rotate=True) plt.gcf().suptitle(f"{p2.name} (coarse alignment)", fontsize=8) plt.gca().set_title(f"{p1.name} - {p2.name}", fontsize=6) - + if viz_coarse_napari: _ = viz_coarse(r1, r2, LEVEL, LEVEL, channel1, channel2, aligner.affine_matrix) if not only_coarse: aligner.ref_img = r1.read_level_channels(LEVEL, channel1) aligner.moving_img = r2.read_level_channels(LEVEL, channel2) - + aligner.compute_shifts() - + fig = aligner.plot_shifts() fig.suptitle(f"{p2.name} (block shift distance)", fontsize=8) fig.axes[0].set_title(p1.name, fontsize=6) @@ -60,9 +60,9 @@ def align_he( set_matplotlib_font(font_size=8) save_all_figs(out_dir=out_dir / 'qc', format='png') - + if not only_qc: - mx = aligner.affine_matrix + mx = aligner.affine_matrix if not only_coarse: mx = aligner.block_affine_matrices_da mosaic = palom.align.block_affine_transformed_moving_img( @@ -84,7 +84,7 @@ def viz_coarse(r1, r2, level1, level2, channel1, channel2, mx): try: import napari except ImportError: - return + return import dask.array as da v = napari.Viewer() is_bf1 = palom.img_util.is_brightfield_img(r1.pyramid[-1][channel1]) @@ -139,7 +139,7 @@ def save_all_figs(dpi=300, format='pdf', out_dir=None): try: import napari except ImportError: print("napari is not installed") else: napari.run() - + ''' Example 1: inspect coarse alignment using napari python align_he.py \ @@ -154,6 +154,6 @@ def save_all_figs(dpi=300, format='pdf', out_dir=None): "X:\crc-scans\histowiz scans\20230105-orion_2_cycles\22199$P54_33_HE$US$SCAN$OR$001 _104050.svs" \ "X:\crc-scans\histowiz scans\20230105-orion_2_cycles\test" \ --px_size1 0.325 - + ''' \ No newline at end of file diff --git a/palom/cli/flow.py b/palom/cli/flow.py index ac766ac..c4bc208 100644 --- a/palom/cli/flow.py +++ b/palom/cli/flow.py @@ -61,7 +61,7 @@ def optical_flow( if mask is None: mask = np.ones(mask_shape, dtype=bool) assert mask.shape == mask_shape - + def func( p1, p2, ch1, ch2, mask, block_info=None @@ -74,7 +74,7 @@ def func( slice(*np.ceil(np.divide(cloc, block_size)).astype(int)) ] return block_optical_flow(chunk1, chunk2, block_size, mask=chunk_mask) - + shifts = da.map_blocks( func, p1=p1, p2=p2, ch1=ch1, ch2=ch2, mask=mask, @@ -104,12 +104,12 @@ def block_optical_flow( wv_img1 = skimage.util.view_as_windows(img1, block_shape, block_shape) wv_img2 = skimage.util.view_as_windows(img2, block_shape, block_shape) h, w = wv_img1.shape[:2] - + wv_mask = np.ones((h, w, 1, 1), dtype=bool) if mask is not None: wv_mask = mask[..., np.newaxis, np.newaxis] wv_mask = wv_mask[..., :h, :w] - + out = np.zeros((3, h, w), dtype=np.float32) product_func = itertools.product if show_progress: @@ -158,7 +158,7 @@ def compose_thumbnail(p1, p2, ch1, ch2, log_intensity=True): rimg2 = skimage.exposure.match_histograms(rimg2, rimg1) out = np.array([rimg1, rimg2, rimg1]) return out - + def plot_legend( shifts, max_radius, plot_scatter, plot_kde, @@ -252,7 +252,7 @@ def process_img_channel_pair( mask = None if compute_mask: mask = reader_block_mask(reader, block_size) - + num_cpus = dask.system.cpu_count() if num_workers > num_cpus: num_workers =num_cpus @@ -309,7 +309,7 @@ def process_img_channel_pair( for ax in fig.get_axes(): ax.set_anchor('N') - + plt.tight_layout() out_dir = img_path.parent / 'qc' out_dir.mkdir(exist_ok=True) @@ -318,7 +318,7 @@ def process_img_channel_pair( bbox_inches='tight', dpi=144 ) figures.append(fig) - + if as_script: return return all_shifts, figures @@ -354,15 +354,15 @@ def process_img_channel_pair( -a, --as_script=AS_SCRIPT Default: True - --- - + --- + NOTE: OTHER_CHANNELS is a list of integers, use [1] to indicate the second channel and [1,2,3] (do not include space) for second, third and fourth channels. Ref https://google.github.io/python-fire/guide/#argument-parsing - - Examples + + Examples python flow.py \ "Z:\RareCyte-S3\P54_CRCstudy_Bridge\S32-Tonsil-P54_Strip_P76.ome.tif" \ diff --git a/palom/cli/helper.py b/palom/cli/helper.py index add56cc..9baf315 100644 --- a/palom/cli/helper.py +++ b/palom/cli/helper.py @@ -11,7 +11,7 @@ def main(argv=sys.argv): 'Generate configuration yaml file for palom-svs run' ) ) - + parser.add_argument( '-i', metavar='input-dir', @@ -36,9 +36,9 @@ def main(argv=sys.argv): help='full path to the resulting configuration yaml file', required=True ) - + args = parser.parse_args(argv[1:]) - + if len(argv) == 1: parser.print_help() return 0 @@ -68,7 +68,7 @@ def main(argv=sys.argv): ) channel_names = [ - '-'.join(p.name.split('_')[-2:][::-1]).replace('.svs', '') + '-'.join(p.name.split('_')[-2:][::-1]).replace('.svs', '') for p in ([ref_slide] + svs_paths) ] diff --git a/palom/cli/schema/svs-config-schema.yml b/palom/cli/schema/svs-config-schema.yml index 7317c48..18c0fb7 100644 --- a/palom/cli/schema/svs-config-schema.yml +++ b/palom/cli/schema/svs-config-schema.yml @@ -31,7 +31,7 @@ image collection: any(include('bright-field image'), include('fluorescent image' _image_base: &_image_base # Filename or the remaining path to the imaging data relative to the input - # dir field + # dir field filename: str() # [optional but recommended] # Channel name in the output ome-tiff @@ -43,7 +43,7 @@ bright-field image: channel names: list(str(), min=1, max=3, required=False) # Output color mode, currently the following 4 modes are available output mode: enum('aec', 'hematoxylin', 'color', 'grayscale', 'dab') - + fluorescent image: <<: *_image_base output mode: str(equals='multichannel') diff --git a/palom/cli/svs.py b/palom/cli/svs.py index b3924cb..6dce0ae 100644 --- a/palom/cli/svs.py +++ b/palom/cli/svs.py @@ -67,11 +67,11 @@ def main(argv=sys.argv): parser_run.set_defaults(func=run) args = parser.parse_args(argv[1:]) - + if len(argv) == 1: parser.print_help() return 0 - + if args.version: print(f"palom v{_version}") return 0 @@ -110,7 +110,7 @@ def run(args): config_data = yamale.make_data(config_file.name) if validate_config(config_data) == 1: return 1 - + config = config_data[0][0] LEVEL = 0 @@ -127,7 +127,7 @@ def run(args): DOWNSCALE_FACTOR = config['pyramid downscale factor'] images = get_image_list(config) - + image_paths = [ pathlib.Path(config['input dir']) / pathlib.Path(i['filename']) for i in images @@ -166,7 +166,7 @@ def run(args): downscale_factor=DOWNSCALE_FACTOR ) - logger.info(f"Finishing {config_file.name}") + logger.info(f"Finishing {config_file.name}") logger.remove(logger_id) # Can't `pathlib.Path.replace` across different disk drive shutil.copy(log_path, qc_path / f"{output_path.name}.log") @@ -209,8 +209,8 @@ def run_palom( ) aligner.coarse_register_affine() - - # FIXME move the saving figure logic + + # FIXME move the saving figure logic plt.suptitle(f"L: {ref_reader.path.name}\nR: {p.name}") fig_w = max(plt.gca().get_xlim()) fig_h = max(plt.gca().get_ylim()) + 100 @@ -219,12 +219,12 @@ def run_palom( plt.tight_layout() plt.savefig(qc_path / f"{idx+1:02d}-{p.name}.png", dpi=72) plt.close() - + aligner.compute_shifts() aligner.constrain_shifts() block_affines.append(aligner.block_affine_matrices_da) - + mosaics = [] m_ref = ref_color_proc.get_processed_color(level=level, mode=img_modes[0]) mosaics.append(m_ref) diff --git a/palom/color.py b/palom/color.py index 683a298..1afe9a8 100644 --- a/palom/color.py +++ b/palom/color.py @@ -22,13 +22,13 @@ def __init__( def find_processed_color_contrast_range(self, mode: str): supported_modes = ['grayscale', 'hematoxylin', 'aec', 'dab'] assert mode in supported_modes - + mode_range = f'_{mode}_range' if hasattr(self, mode_range): return getattr(self, mode_range) - + contrast_ref_img = self.contrast_img - + process_func = self.__getattribute__(f"rgb2{mode}") img = process_func(contrast_ref_img) @@ -68,7 +68,7 @@ def get_processed_color(self, rgb_img, mode='grayscale'): process_func = self.__getattribute__(f"rgb2{mode}") intensity_range = self.find_processed_color_contrast_range(mode) - + processed = da.map_blocks( process_func, rgb_img, @@ -96,7 +96,7 @@ def __init__( thumbnail_level = len(pyramid) - 1 super().__init__(pyramid[thumbnail_level], channel_axis=0) self.pyramid = pyramid - + def get_processed_color(self, level, mode='grayscale', out_dtype=None): if mode == 'color': return self.pyramid[level] diff --git a/palom/extract.py b/palom/extract.py index 5d241e8..4696fd9 100644 --- a/palom/extract.py +++ b/palom/extract.py @@ -15,7 +15,7 @@ def imagej_rgb2cmyk(rgb_img): cmy = 1-rgb k = cmy.min(axis=0) - + s = 1-k with warnings.catch_warnings(): diff --git a/palom/img_util.py b/palom/img_util.py index 4917147..6beac37 100644 --- a/palom/img_util.py +++ b/palom/img_util.py @@ -18,7 +18,7 @@ def cv2_pyramid(img, max_size=1024): def whiten(img, sigma=1): border_mode = cv2.BORDER_REFLECT g_img = img if sigma == 0 else cv2.GaussianBlur( - img, (0, 0), sigma, + img, (0, 0), sigma, borderType=border_mode ).astype(np.float32) log_img = cv2.Laplacian( @@ -80,25 +80,25 @@ def block_labeled_mask(img_shape, block_shape, out_chunks=None): return da.from_array(full_mask, chunks=out_chunks) -def to_napari_affine(mx): +def to_napari_affine(mx): ul = np.flip(mx[:2, :2], (0, 1)) rows = np.hstack([ul, np.flipud(mx[:2, 2:3])]) return np.vstack((rows, [0, 0, 1])) # orders of magnitute faster than skimage.transform.downscale_local_mean -# also the gives sensible values of pixels on the edge -def cv2_downscale_local_mean(img, factor): - assert img.ndim in [2, 3] - img = np.asarray(img) - axis_moved = False - channel_ax = np.argmin(img.shape) - if (img.ndim == 3) & (channel_ax != 2): - img = np.moveaxis(img, channel_ax, 2) +# also the gives sensible values of pixels on the edge +def cv2_downscale_local_mean(img, factor): + assert img.ndim in [2, 3] + img = np.asarray(img) + axis_moved = False + channel_ax = np.argmin(img.shape) + if (img.ndim == 3) & (channel_ax != 2): + img = np.moveaxis(img, channel_ax, 2) axis_moved = True simg = cv2.blur( img, ksize=(factor, factor), anchor=(0, 0) )[::factor, ::factor] - if axis_moved: - simg = np.moveaxis(simg, 2, channel_ax) + if axis_moved: + simg = np.moveaxis(simg, 2, channel_ax) return simg \ No newline at end of file diff --git a/palom/pyramid.py b/palom/pyramid.py index ac0a1ea..faf3e96 100644 --- a/palom/pyramid.py +++ b/palom/pyramid.py @@ -135,10 +135,10 @@ def write_pyramid( tile_shapes = pyramid_setting.tile_shapes(base_shape) shapes = pyramid_setting.pyramid_shapes(base_shape) - if tile_size is not None: - assert tile_size % 16 == 0, ( - f"tile_size must be None or multiples of 16, not {tile_size}" - ) + if tile_size is not None: + assert tile_size % 16 == 0, ( + f"tile_size must be None or multiples of 16, not {tile_size}" + ) tile_shapes = [(tile_size, tile_size)] * num_levels dtype = ref_m.dtype @@ -214,7 +214,7 @@ def count_num_channels(imgs): 1 if img.ndim == 2 else img.shape[0] for img in imgs ]) - + def tile_from_combined_mosaics(mosaics, tile_shape, save_RAM=False): num_rows, num_cols = mosaics[0].shape[1:3] @@ -259,7 +259,7 @@ def tile_from_pyramid( # img = tifffile.imread(path, series=0, level=level, key=c) if not is_mask: img = img.map_blocks( - cv2.blur, + cv2.blur, ksize=(downscale_factor, downscale_factor), anchor=(0, 0) ) img = img.persist() if save_RAM else img.compute() diff --git a/palom/reader.py b/palom/reader.py index 0ffc745..9a7804f 100644 --- a/palom/reader.py +++ b/palom/reader.py @@ -30,11 +30,11 @@ def validate_pyramid(pyramid: list[da.Array], channel_axis:int) -> bool: assert level.ndim == 3, '' if np.argmin(level.shape) != channel_axis: logger.warning( - f"level {i} has shape of {level.shape} while given" + f"level {i} has shape of {level.shape} while given" f" `channel_axis` is {channel_axis}" ) return True - + def normalize_axis_order(self): if self.channel_axis == 0: return self.pyramid @@ -83,7 +83,7 @@ def level_downsamples(self) -> dict[int, int]: i: round(self.pyramid[0].shape[1] / level.shape[1]) for i, level in enumerate(self.pyramid) } - + @property def pixel_dtype(self) -> np.dtype: return self.pyramid[0].dtype @@ -151,7 +151,7 @@ def __init__( # https://python-poetry.org/docs/pyproject/#extras # https://github.com/AllenCellModeling/aicsimageio/blob/main/aicsimageio/readers/bioformats_reader.py#L33-L40 from napari_lazy_openslide import OpenSlideStore - + self.path = pathlib.Path(path) self.store = OpenSlideStore(str(self.path)) self.zarr = zarr.open(self.store, mode='r') @@ -165,7 +165,7 @@ def pyramid_from_svs(self) -> list[da.Array]: da.from_zarr(self.store, component=d['path'])[..., :3] for d in self.zarr.attrs['multiscales'][0]['datasets'] ] - + @property def pixel_size(self): if self._pixel_size is not None: diff --git a/palom/register.py b/palom/register.py index bb8178e..6491e59 100644 --- a/palom/register.py +++ b/palom/register.py @@ -19,15 +19,15 @@ register_translation = skimage.feature.register_translation -# +# # Image-based registration -# +# def phase_cross_correlation(img1, img2, sigma, upsample=10, module='skimage'): assert module in ['cv2', 'skimage'] - + img1w = img_util.whiten(img1, sigma) img2w = img_util.whiten(img2, sigma) - + if module == 'skimage': with warnings.catch_warnings(): warnings.filterwarnings( @@ -44,7 +44,7 @@ def phase_cross_correlation(img1, img2, sigma, upsample=10, module='skimage'): shift, _error, _phasediff = register_translation( img1w, img2w, **kwargs ) - + elif module == 'cv2': shift_xy, _response = cv2.phaseCorrelate(img1w, img2w) shift = shift_xy[::-1] @@ -75,7 +75,7 @@ def cv2_translate(img, shift): assert img.ndim == len(shift) == 2 sy, sx = shift return cv2.warpAffine( - img, + img, np.array([[1, 0, sx], [0, 1, sy]], dtype=float), img.shape[::-1] ) @@ -91,9 +91,9 @@ def normalized_phase_correlation(img1, img2, sigma): return corr -# +# # Feature-based registration -# +# def feature_based_registration( img_left, img_right, detect_flip_rotate=False, @@ -103,7 +103,7 @@ def feature_based_registration( flip_rotate_func, mx_fr = np.array, np.eye(3) if detect_flip_rotate: flip_rotate_func, mx_fr = match_test_flip_rotate(img_left, img_right) - + img_right = flip_rotate_func(img_right) mx_affine = ensambled_match( img_left, img_right, @@ -180,7 +180,7 @@ def ensambled_match( all_dst = np.vstack([i[1] for i in all_found]) t_matrix, mask = cv2.estimateAffine2D( - all_dst, all_src, + all_dst, all_src, method=cv2.RANSAC, ransacReprojThreshold=ransacReprojThreshold, maxIters=5000 @@ -231,7 +231,7 @@ def cv2_feature_detect_and_match( [keypoints_right[m.trainIdx].pt for m in matches] ) t_matrix, mask = cv2.estimateAffine2D( - dst_pts, src_pts, + dst_pts, src_pts, method=cv2.RANSAC, ransacReprojThreshold=30, maxIters=5000 ) if plot_match_result == True: diff --git a/palom/register_util.py b/palom/register_util.py index 0c18f10..871bd8f 100644 --- a/palom/register_util.py +++ b/palom/register_util.py @@ -42,7 +42,7 @@ def match_bf_fl_histogram(img1, img2): # TODO does it make a difference to min/max rescale before histogram # matching? is_bf_img1, is_bf_img2 = [ - img_util.is_brightfield_img(i) + img_util.is_brightfield_img(i) for i in (img1, img2) ] if is_bf_img1 == is_bf_img2: @@ -61,7 +61,7 @@ def plot_img_keypoints(imgs, keypoints): flags=cv2.DRAW_MATCHES_FLAGS_DEFAULT )) a.set_title(len(k)) - return + return def get_flip_mx(img_shape, flip_axis):