diff --git a/src/classpose/entrypoints/predict_wsi.py b/src/classpose/entrypoints/predict_wsi.py index bb57a8d..23bf8d6 100644 --- a/src/classpose/entrypoints/predict_wsi.py +++ b/src/classpose/entrypoints/predict_wsi.py @@ -146,6 +146,8 @@ def __init__( self.ts = manager.Value("f", 0) self.mpp_x = manager.Value("f", 0) self.mpp_y = manager.Value("f", 0) + self.bounds_x = manager.Value("f", 0) + self.bounds_y = manager.Value("f", 0) self.tissue_cnts = manager.list() self.roi_cnts = manager.list() @@ -169,6 +171,17 @@ def _init_slide(self): self.mpp = get_slide_resolution(self.slide) self.mpp_x.value = self.mpp[0] self.mpp_y.value = self.mpp[1] + bounds_x = self.slide.properties.get("openslide.bounds-x") + bounds_y = self.slide.properties.get("openslide.bounds-y") + self.bounds_x.value = float(bounds_x) if bounds_x is not None else 0.0 + self.bounds_y.value = float(bounds_y) if bounds_y is not None else 0.0 + logger.info(f"Slide bounds offset: x={self.bounds_x.value}, y={self.bounds_y.value}") + + if self.roi_tree is not None and ( + self.bounds_x.value != 0 or self.bounds_y.value != 0 + ): + self._align_roi_tree_to_slide_bounds() + target_downsample = min( self.train_mpp / self.mpp[0], self.train_mpp / self.mpp[1] ) @@ -196,6 +209,27 @@ def _init_slide(self): logger.info(f"Slide downsample: {self.ts.value}") logger.info(f"Level: {self.level}") + def _align_roi_tree_to_slide_bounds(self): + """ + Align ROI polygons to slide coordinates when OpenSlide bounds offsets are present. + + QuPath ROI annotations are exported in coordinates relative to the displayed image + origin (i.e. bounds-offset corrected). Internal OpenSlide tile coordinates are in + level-0 slide space. So if needed, shift ROI polygons by +bounds offsets so + ROI-based tile selection/filtering uses the same frame. + """ + bounds_x_val = float(self.bounds_x.value) + bounds_y_val = float(self.bounds_y.value) + geometries = list(self.roi_tree.geometries) + + logger.info(f"Applying bounds offset to ROI polygons: x={bounds_x_val}, y={bounds_y_val}") + + shifted = [ + shapely.affinity.translate(geom, xoff=bounds_x_val, yoff=bounds_y_val) + for geom in geometries + ] + self.roi_tree = shapely.STRtree(shifted) + def _get_tissue_contours(self): if self.tissue_detection_model_path is not None: logger.info("Detecting tissue contours using GrandQC") @@ -460,7 +494,7 @@ def __call__( continue center = np.round(polygon.centroid.coords[0], 2).tolist() curr_coords = curr_coords.tolist() - curr_coords.append(curr_coords[0]) + curr_coords.append(curr_coords[0].copy()) cl = class_masks[cell_mask][0] curr_cell = { "id": str(uuid.uuid4()), @@ -611,6 +645,44 @@ def to_geojson_polygon(curr_cell: dict) -> dict: return curr_cell +def apply_bounds_offset_to_feature(feature: dict, bounds_x: float, bounds_y: float) -> dict: + """ + Apply bounds offset to a GeoJSON Polygon feature's geometry and centroid measurements. + + Args: + feature: GeoJSON Feature dict with Polygon geometry + bounds_x: X offset from openslide.bounds-x property + bounds_y: Y offset from openslide.bounds-y property + + Returns: + Feature with shifted polygon coordinates and updated centroid measurements. + """ + if not feature or "geometry" not in feature: + return feature + + geometry = feature["geometry"] + if "coordinates" not in geometry: + return feature + + shifted_coordinates = [] + for ring in geometry["coordinates"]: + shifted_ring = [ + [point[0] - bounds_x, point[1] - bounds_y] + for point in ring + ] + shifted_coordinates.append(shifted_ring) + geometry["coordinates"] = shifted_coordinates + + if "properties" in feature and "measurements" in feature["properties"]: + for measurement in feature["properties"]["measurements"]: + if measurement["name"] == "centroidX": + measurement["value"] -= bounds_x + elif measurement["name"] == "centroidY": + measurement["value"] -= bounds_y + + return feature + + def deduplicate(features: list[dict], max_dist: float = 15 / 2) -> list[dict]: """ Deduplicate cells based on their size and location. @@ -705,7 +777,7 @@ def shapely_polygon_to_geojson( coords = np.array(polygon.exterior.coords.xy).T center = coords.mean(axis=1).tolist() coords = coords.tolist() - coords.append(coords[0]) + coords.append(coords[0].copy()) feature_id = str(uuid.uuid4()) properties = { "objectType": object_type, @@ -1241,6 +1313,15 @@ def main(args): logger.info("Filtering cells based on tissue contours") tissue_cnts = list(slide.tissue_cnts) polygons = filter_cells_by_contours(polygons, tissue_cnts) + + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + tissue_cnts = [ + shapely.affinity.translate(cnt, xoff=-bounds_x_val, yoff=-bounds_y_val) + for cnt in tissue_cnts + ] + tissue_features = [] for i, cnt in enumerate(tissue_cnts): tissue_features.extend( @@ -1311,6 +1392,15 @@ def main(args): if polygon is not None: artefact_polygons.append(polygon) + + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + artefact_polygons = [ + shapely.affinity.translate(poly, xoff=-bounds_x_val, yoff=-bounds_y_val) + for poly in artefact_polygons + ] + artefact_features = [] for i, poly in enumerate(artefact_polygons): artefact_features.extend( @@ -1340,6 +1430,15 @@ def main(args): with open(output_folder / artefact_contours_filename, "w") as f: json.dump(artefact_contours_fmt, f) + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + logger.info(f"Applying bounds offset to output coordinates: x={bounds_x_val}, y={bounds_y_val}") + polygons = [ + apply_bounds_offset_to_feature(f, bounds_x_val, bounds_y_val) + for f in polygons + ] + geojson_fmt = { "type": "FeatureCollection", "features": polygons, diff --git a/src/classpose/entrypoints/predict_wsi_cpsam.py b/src/classpose/entrypoints/predict_wsi_cpsam.py index cee161c..70cc52b 100644 --- a/src/classpose/entrypoints/predict_wsi_cpsam.py +++ b/src/classpose/entrypoints/predict_wsi_cpsam.py @@ -72,6 +72,7 @@ DEFAULT_TRAIN_MPP, MIN_TILE_SIZE, SlideLoader, + apply_bounds_offset_to_feature, deduplicate, filter_cells_by_artefacts, filter_cells_by_contours, @@ -434,6 +435,17 @@ def main(args): logger.info("Filtering cells based on tissue contours") tissue_cnts = list(slide.tissue_cnts) polygons = filter_cells_by_contours(polygons, tissue_cnts) + + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + tissue_cnts = [ + shapely.affinity.translate( + cnt, xoff=-bounds_x_val, yoff=-bounds_y_val + ) + for cnt in tissue_cnts + ] + tissue_features = [] for i, cnt in enumerate(tissue_cnts): tissue_features.extend( @@ -476,7 +488,7 @@ def main(args): artefact_cnts, artefact_geojson, ) = detect_artefacts_wsi( - slide=OpenSlide(args.slide_path), + slide=OpenSlide(slide.real_slide_path), model_art_path=args.artefact_detection_model_path, model_td_path=args.tissue_detection_model_path, device=devices[0], @@ -508,6 +520,16 @@ def main(args): if polygon is not None: artefact_polygons.append(polygon) + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + artefact_polygons = [ + shapely.affinity.translate( + poly, xoff=-bounds_x_val, yoff=-bounds_y_val + ) + for poly in artefact_polygons + ] + artefact_features = [] for i, poly in enumerate(artefact_polygons): artefact_features.extend( @@ -539,6 +561,19 @@ def main(args): with open(output_folder / artefact_contours_filename, "w") as f: json.dump(artefact_contours_fmt, f) + bounds_x_val = float(slide.bounds_x.value) + bounds_y_val = float(slide.bounds_y.value) + if bounds_x_val != 0 or bounds_y_val != 0: + logger.info( + "Applying bounds offset to output coordinates: x=%s, y=%s", + bounds_x_val, + bounds_y_val, + ) + polygons = [ + apply_bounds_offset_to_feature(f, bounds_x_val, bounds_y_val) + for f in polygons + ] + geojson_fmt = { "type": "FeatureCollection", "features": polygons, diff --git a/src/classpose/grandqc/wsi_artefact_detection.py b/src/classpose/grandqc/wsi_artefact_detection.py index 7b02074..37707bf 100644 --- a/src/classpose/grandqc/wsi_artefact_detection.py +++ b/src/classpose/grandqc/wsi_artefact_detection.py @@ -62,6 +62,7 @@ def detect_artefacts_wsi( mpp_model_td: int = 10, m_p_s_model_td: int = 512, min_area: int = 0, + apply_bounds_offset: bool = False, ) -> tuple[np.ndarray, np.ndarray, dict[str, Any], dict[str, Any]]: """ Detects artefacts in a whole-slide image using GrandQC method. @@ -85,6 +86,9 @@ def detect_artefacts_wsi( mpp_model_td (int): MPP for tissue detection. m_p_s_model_td (int): patch size for tissue detection. min_area (int): minimum area for tissue polygons. + apply_bounds_offset (bool): if True, shift artefact contours and + GeoJSON by OpenSlide bounds offsets, so output coordinates are + relative to the displayed image origin. Defaults to False. Returns: tuple: artefact_mask (np.ndarray), artefact_map (np.ndarray), artefact_cnts (dict), geojson (dict) @@ -100,8 +104,12 @@ def detect_artefacts_wsi( encoder_model_weights=encoder_model_weights, device=device, min_area=min_area, + apply_bounds_offset=False, ) + bounds_x = float(slide.properties.get("openslide.bounds-x", 0.0)) + bounds_y = float(slide.properties.get("openslide.bounds-y", 0.0)) + # Load artefact model model_art_path = download_if_unavailable(model_art_path, MODEL_URL_PATH) preprocessing_fn = smp.encoders.get_preprocessing_fn( @@ -306,6 +314,27 @@ def detect_artefacts_wsi( min_artifact_area, ) + if apply_bounds_offset and (bounds_x != 0 or bounds_y != 0): + grandqc_logger.info( + "Applying bounds offset to artefact output coordinates: x=%s, y=%s", + bounds_x, + bounds_y, + ) + + for cnt in artefact_cnts.values(): + cnt["contour"] = cnt["contour"] - np.array([bounds_x, bounds_y]) + if cnt["holes"]: + cnt["holes"] = [ + hole - np.array([bounds_x, bounds_y]) for hole in cnt["holes"] + ] + + for feature in geojson["features"]: + rings = feature["geometry"]["coordinates"] + feature["geometry"]["coordinates"] = [ + [[point[0] - bounds_x, point[1] - bounds_y] for point in ring] + for ring in rings + ] + del model del preprocessing_fn del slide @@ -370,6 +399,7 @@ def detect_artefacts_wsi( model_td_path=args.model_td_path, min_area=args.min_area, device=args.device, + apply_bounds_offset=True, ) # Save outputs diff --git a/src/classpose/grandqc/wsi_tissue_detection.py b/src/classpose/grandqc/wsi_tissue_detection.py index 7dec55e..926900b 100644 --- a/src/classpose/grandqc/wsi_tissue_detection.py +++ b/src/classpose/grandqc/wsi_tissue_detection.py @@ -35,6 +35,7 @@ def detect_tissue_wsi( device: str = "cuda", # filtering parameters min_area: int = 0, + apply_bounds_offset: bool = False, ) -> tuple[Image, np.ndarray, np.ndarray, list[np.ndarray], dict[str, Any],]: """ Detects tissue in a whole-slide image using a U-Net model. This is adapted @@ -58,6 +59,9 @@ def detect_tissue_wsi( device (str, optional): device to run the model on. Defaults to "cuda". min_area (int, optional): minimum area of the tissue polygons. Defaults to 0. + apply_bounds_offset (bool, optional): if True, shift output contours and + GeoJSON by OpenSlide bounds offsets, so output coordinates are relative + to the displayed image origin. Defaults to False. Returns: tuple: tuple containing the following elements: @@ -86,6 +90,9 @@ def detect_tissue_wsi( model.to(device) model.eval() + bounds_x = float(slide.properties.get("openslide.bounds-x", 0.0)) + bounds_y = float(slide.properties.get("openslide.bounds-y", 0.0)) + w_l0, h_l0, mpp, thumbnail_dimensions = extract_slide_info( slide, mpp_model_td ) @@ -265,6 +272,26 @@ def detect_tissue_wsi( } ) + if apply_bounds_offset and (bounds_x != 0 or bounds_y != 0): + grandqc_logger.info( + "Applying bounds offset to tissue output coordinates: x=%s, y=%s", + bounds_x, + bounds_y, + ) + for cnt in output_cnts.values(): + cnt["contour"] = cnt["contour"] - np.array([bounds_x, bounds_y]) + if cnt["holes"]: + cnt["holes"] = [ + hole - np.array([bounds_x, bounds_y]) for hole in cnt["holes"] + ] + + for feature in geojson["features"]: + rings = feature["geometry"]["coordinates"] + feature["geometry"]["coordinates"] = [ + [[point[0] - bounds_x, point[1] - bounds_y] for point in ring] + for ring in rings + ] + del model del preprocessing_fn del slide @@ -322,6 +349,7 @@ def detect_tissue_wsi( model_td_path=args.model_path, min_area=args.min_area, device=args.device, + apply_bounds_offset=True, ) image.save(args.output_path + "_image.png")