diff --git a/CHANGELOG.md b/CHANGELOG.md index 4a576f7..edbfb35 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,2 +1,8 @@ +This file might no be up to date nor complete. Please check the Releases page for more information on versions. + +## Version 2024.12.10 +- Fixed offset writing geojson file in segmentation module. +- Added original_pixelsize parameters to rescale coordinates to match final image size. This allows the use of QuPath pixel classifier trained on resized image, eg. with a lower Resolution parameter (higher pixel size). + ## Version 2024.11.19 - Initial public release. \ No newline at end of file diff --git a/histoquant/seg.py b/histoquant/seg.py index 548d166..5882b4e 100644 --- a/histoquant/seg.py +++ b/histoquant/seg.py @@ -32,7 +32,11 @@ def get_pixelsize(image_name: str) -> float: """ with tifffile.TiffFile(image_name) as tif: - return 1e6 / tif.pages[0].tags["XResolution"].value[0] + # XResolution is a tuple, numerator, denomitor. The inverse is the pixel size + return ( + tif.pages[0].tags["XResolution"].value[1] + / tif.pages[0].tags["XResolution"].value[0] + ) def convert_to_pixels(filters, pixelsize): @@ -142,7 +146,7 @@ def get_image_skeleton(img: np.ndarray, minsize=0) -> np.ndarray: def get_collection_from_skel( - skeleton: Skeleton, properties: dict + skeleton: Skeleton, properties: dict, rescale_factor: float = 1.0, offset=0.5 ) -> geojson.FeatureCollection: """ Get the coordinates of each skeleton path as a GeoJSON Features in a @@ -154,6 +158,11 @@ def get_collection_from_skel( skeleton : skan.Skeleton properties : dict QuPatj objects' properties. + rescale_factor : float + Rescale output coordinates by this factor. + offset : float + Shift coordinates by this amount, typically to get pixel centers or edges. + Default is 0.5. Returns ------- @@ -171,7 +180,7 @@ def get_collection_from_skel( collection.append( geojson.Feature( geometry=shapely.LineString( - skeleton.path_coordinates(ind)[:, ::-1] + (skeleton.path_coordinates(ind)[:, ::-1] + offset) * rescale_factor ), # shape object properties=prop, # object properties id=str(uuid.uuid4()), # object uuid @@ -182,7 +191,7 @@ def get_collection_from_skel( def get_collection_from_poly( - contours: list, properties: dict + contours: list, properties: dict, rescale_factor: float = 1.0, offset: float = 0.5 ) -> geojson.FeatureCollection: """ Gather coordinates in the list and put them in GeoJSON format as Polygons. @@ -195,6 +204,11 @@ def get_collection_from_poly( contours : list properties : dict QuPatj objects' properties. + rescale_factor : float + Rescale output coordinates by this factor. + offset : float + Shift coordinates by this amount, typically to get pixel centers or edges. + Default is 0.5. Returns ------- @@ -202,10 +216,11 @@ def get_collection_from_poly( A FeatureCollection ready to be written as geojson. """ - collection = [ geojson.Feature( - geometry=shapely.Polygon(np.fliplr(contour)), # shape object + geometry=shapely.Polygon( + np.fliplr((contour + offset) * rescale_factor) + ), # shape object properties=properties, # object properties id=str(uuid.uuid4()), # object uuid ) @@ -216,7 +231,7 @@ def get_collection_from_poly( def get_collection_from_points( - coords: list, properties: dict + coords: list, properties: dict, rescale_factor: float = 1.0, offset: float = 0.5 ) -> geojson.FeatureCollection: """ Gather coordinates from `coords` and put them in GeoJSON format. @@ -228,6 +243,8 @@ def get_collection_from_points( ---------- coords : list properties : dict + rescale_factor : float + Rescale output coordinates by this factor. Returns ------- @@ -237,7 +254,9 @@ def get_collection_from_points( collection = [ geojson.Feature( - geometry=shapely.Point(np.flip(coord)), # shape object + geometry=shapely.Point( + np.flip((coord + offset) * rescale_factor) + ), # shape object properties=properties, # object properties id=str(uuid.uuid4()), # object uuid ) @@ -248,7 +267,7 @@ def get_collection_from_points( def segment_lines( - img: np.ndarray, geojson_props: dict, minsize=0.0 + img: np.ndarray, geojson_props: dict, minsize=0.0, rescale_factor=1.0 ) -> geojson.FeatureCollection: """ Wraps skeleton analysis to get paths coordinates. @@ -261,6 +280,8 @@ def segment_lines( GeoJSON properties of objects. minsize : float Minimum size in pixels for an object. + rescale_factor : float + Rescale output coordinates by this factor. Returns ------- @@ -273,7 +294,9 @@ def segment_lines( # get paths coordinates as FeatureCollection skeleton = Skeleton(skel, keep_images=False) - return get_collection_from_skel(skeleton, geojson_props) + return get_collection_from_skel( + skeleton, geojson_props, rescale_factor=rescale_factor + ) def segment_polygons( @@ -283,6 +306,7 @@ def segment_polygons( area_max: float = np.inf, ecc_min: float = 0.0, ecc_max: float = 1.0, + rescale_factor: float = 1.0, ) -> geojson.FeatureCollection: """ Polygon segmentation. @@ -297,6 +321,8 @@ def segment_polygons( Minimum and maximum area in pixels for an object. ecc_min, ecc_max : float Minimum and maximum eccentricity for an object. + rescale_factor: float + Rescale output coordinates by this factor. Returns ------- @@ -328,7 +354,9 @@ def segment_polygons( label_image = label_image > 0 contours = measure.find_contours(label_image) - return get_collection_from_poly(contours, geojson_props) + return get_collection_from_poly( + contours, geojson_props, rescale_factor=rescale_factor + ) def segment_points( @@ -339,6 +367,7 @@ def segment_points( ecc_min: float = 0, ecc_max: float = 1, dist_thresh: float = 0, + rescale_factor: float = 1, ) -> geojson.FeatureCollection: """ Point segmentation. @@ -359,6 +388,8 @@ def segment_points( dist_thresh : float Maximal distance in pixels between objects before considering them as isolated and remove them. 0 disables it. + rescale_factor : float + Rescale output coordinates by this factor. Returns ------- @@ -410,4 +441,6 @@ def segment_points( # get points coordinates coords = np.argwhere(bw) - return get_collection_from_points(coords, geojson_props) + return get_collection_from_points( + coords, geojson_props, rescale_factor=rescale_factor + ) diff --git a/pyproject.toml b/pyproject.toml index a720b05..66ad841 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "histoquant" -version = "2024.11.27" +version = "2024.12.10" authors = [{ name = "Guillaume Le Goc", email = "g.legoc@posteo.org" }] description = "Quantification of objects in histological slices" readme = "README.md" diff --git a/scripts/segmentation/segment_images.py b/scripts/segmentation/segment_images.py index 5d91e19..8716bdd 100644 --- a/scripts/segmentation/segment_images.py +++ b/scripts/segmentation/segment_images.py @@ -18,8 +18,8 @@ To exclude objects near the edges of an ROI, specify the path to masks stored as images with the same names as probabilities images (without their suffix). -Author : Guillaume Le Goc (g.legoc@posteo.org) @ NeuroPSI -Version : 2024.11.27 +author : Guillaume Le Goc (g.legoc@posteo.org) @ NeuroPSI +version : 2024.12.10 """ @@ -37,17 +37,22 @@ pd.options.mode.copy_on_write = True # prepare for pandas 3 # --- Parameters -IMAGES_DIR = "/path/to/images/to/segment" +IMAGES_DIR = "/path/to/images" """Full path to the images to segment.""" MASKS_DIR = "path/to/corresponding/masks" """Full path to the masks, to exclude objects near the brain edges (set to None or empty string to disable this feature).""" MASKS_EXT = "tiff" """Masks files extension.""" -SEGTYPE = "fibers" +SEGTYPE = "boutons" """Type of segmentation.""" IMG_SUFFIX = "_Probabilities.tiff" """Images suffix, including extension. Masks must be the same name without the suffix.""" +ORIGINAL_PIXELSIZE = 0.4500 +"""Original images pixel size in microns. This is in case the pixel classifier uses +a lower resolution, yielding smaller probability maps, so output objects coordinates +need to be rescaled to the full size images. The pixel size is written in the "Image" +tab in QuPath.""" CHANNELS_PARAMS = [ { @@ -60,7 +65,7 @@ { "name": "dsred", "target_channel": 1, - "proba_threshold": 0.85, + "proba_threshold": 0.65, "qp_class": "Fibers: DsRed", "qp_color": [224, 153, 18], }, @@ -80,13 +85,13 @@ - qp_class: str, name of QuPath classification - qp_color: list of RGB values, associated color""" -EDGE_DIST = 50 +EDGE_DIST = 0 """Distance to brain edge to ignore, in µm. 0 to disable.""" FILTERS = { "length_low": 1.5, # minimal length in microns - for lines - "area_low": 1.1, # minimal area in µm² - for polygons and points - "area_high": 10, # maximal area in µm² - for polygons and points + "area_low": 10, # minimal area in µm² - for polygons and points + "area_high": 1000, # maximal area in µm² - for polygons and points "ecc_low": 0.0, # minimal eccentricity - for polygons and points (0 = circle) "ecc_high": 0.9, # maximal eccentricity - for polygons and points (1 = line) "dist_thresh": 30, # maximal inter-point distance in µm - for points @@ -237,7 +242,9 @@ def parameters_as_dict( } -def write_parameters(outfile: str, parameters: dict, filters: dict): +def write_parameters( + outfile: str, parameters: dict, filters: dict, original_pixelsize: float +): """ Write parameters to `outfile`. @@ -252,12 +259,16 @@ def write_parameters(outfile: str, parameters: dict, filters: dict): General parameters. filters : dict Filters parameters. + original_pixelsize : float + Size of pixels in original image. """ with open(outfile, "w") as fid: fid.writelines(f"date = {datetime.now().strftime('%d-%B-%Y %H:%M:%S')}\n") + fid.writelines(f"original_pixelsize = {original_pixelsize}\n") + for key, value in parameters.items(): fid.writelines(f"{key} = {value}\n") @@ -271,6 +282,7 @@ def process_directory( images_dir: str, img_suffix: str = "", segtype: str = "", + original_pixelsize: float = 1.0, target_channel: int = 0, proba_threshold: float = 0.0, qupath_class: str = "Object", @@ -292,6 +304,8 @@ def process_directory( Images suffix, including extension. segtype : str Segmentation type. + original_pixelsize : float + Original images pixel size in microns. target_channel : int Index of the channel containning the objects of interest (eg. not the background), in the probability map (*not* the original images channels). @@ -339,13 +353,16 @@ def process_directory( if os.path.isfile(param_file): raise FileExistsError("Parameters file already exists.") else: - write_parameters(param_file, parameters, filters) + write_parameters(param_file, parameters, filters, original_pixelsize) - # convert parameters to pixels + # convert parameters to pixels in probability map pixelsize = hq.seg.get_pixelsize(images_list[0]) # get pixel size edge_dist = int(edge_dist / pixelsize) filters = hq.seg.convert_to_pixels(filters, pixelsize) + # get rescaling factor + rescale_factor = pixelsize / original_pixelsize + # get GeoJSON properties geojson_props = get_geojson_properties( qupath_class, qupath_color, objtype=QUPATH_TYPE @@ -385,7 +402,10 @@ def process_directory( if seg_method == "lines": collection = hq.seg.segment_lines( - img, geojson_props, minsize=filters["length_low"] + img, + geojson_props, + minsize=filters["length_low"], + rescale_factor=rescale_factor, ) elif seg_method == "polygons": @@ -396,6 +416,7 @@ def process_directory( area_max=filters["area_high"], ecc_min=filters["ecc_low"], ecc_max=filters["ecc_high"], + rescale_factor=rescale_factor, ) elif seg_method == "points": @@ -407,6 +428,7 @@ def process_directory( ecc_min=filters["ecc_low"], ecc_max=filters["ecc_high"], dist_thresh=filters["dist_thresh"], + rescale_factor=rescale_factor, ) else: # we already printed an error message @@ -445,6 +467,7 @@ def make_suffix(s): IMAGES_DIR, img_suffix=IMG_SUFFIX, segtype=SEGTYPE, + original_pixelsize=ORIGINAL_PIXELSIZE, target_channel=param["target_channel"], proba_threshold=param["proba_threshold"], qupath_class=param["qp_class"],