diff --git a/PyNutil/io/propagation.py b/PyNutil/io/propagation.py index 22a611fe44b7ee85384cc70a56dcc942154ea5cd..600ad6033c5b296b17a6037a183ac15b14d160bb 100644 --- a/PyNutil/io/propagation.py +++ b/PyNutil/io/propagation.py @@ -1,5 +1,6 @@ import math, re + def propagate(arr): arr = arr.copy() for slice in arr: diff --git a/PyNutil/io/reconstruct_dzi.py b/PyNutil/io/reconstruct_dzi.py index 862c2ec55078ed9e0ddbcdfc8f96f48a07ad476b..bbdff4038d0577d2d0fba447d657af634627801f 100644 --- a/PyNutil/io/reconstruct_dzi.py +++ b/PyNutil/io/reconstruct_dzi.py @@ -4,6 +4,7 @@ import os import zipfile import xmltodict + def reconstruct_dzi(zip_file_path): """ Reconstructs a Deep Zoom Image (DZI) from a zip file containing the tiles. diff --git a/PyNutil/main.py b/PyNutil/main.py index 8522cea0c6d5d2e999c3d0c19676ab719d57ee26..3642a72d5a1e98f84e6dfc1828dd4ebe65ec667d 100644 --- a/PyNutil/main.py +++ b/PyNutil/main.py @@ -1,10 +1,15 @@ import json from .io.atlas_loader import load_atlas_data, load_custom_atlas -from .processing.data_analysis import quantify_labeled_points, map_to_custom_regions, apply_custom_regions +from .processing.data_analysis import ( + quantify_labeled_points, + map_to_custom_regions, + apply_custom_regions, +) from .io.file_operations import save_analysis_output from .io.read_and_write import open_custom_region_file from .processing.coordinate_extraction import folder_to_atlas_space + class PyNutil: """ A class used to perform brain-wide quantification and spatial analysis of features in serial section images. @@ -150,7 +155,7 @@ class PyNutil: self.centroids_len, self.segmentation_filenames, self.per_point_undamaged, - self.per_centroid_undamaged + self.per_centroid_undamaged, ) = folder_to_atlas_space( self.segmentation_folder, self.alignment_json, @@ -201,7 +206,7 @@ class PyNutil: self.centroids_labels, self.atlas_labels, self.per_point_undamaged, - self.per_centroid_undamaged + self.per_centroid_undamaged, ) if self.custom_regions_dict is not None: self.custom_label_df, self.label_df = apply_custom_regions( diff --git a/PyNutil/processing/coordinate_extraction.py b/PyNutil/processing/coordinate_extraction.py index ddf569eb565e14624d939792c89b9122ca43d878..f6d26c2bfdd6f8aa7773f6df3f7b2519d3effd63 100644 --- a/PyNutil/processing/coordinate_extraction.py +++ b/PyNutil/processing/coordinate_extraction.py @@ -43,15 +43,17 @@ def get_centroids_and_area(segmentation, pixel_cut_off=0): coords = np.array([label.coords for label in labels_info], dtype=object) return centroids, area, coords + def update_spacing(anchoring, width, height, grid_spacing): if len(anchoring) != 9: print("Anchoring does not have 9 elements.") - ow = np.sqrt(sum([anchoring[i+3] ** 2 for i in range(3)])) - oh = np.sqrt(sum([anchoring[i+6] ** 2 for i in range(3)])) + ow = np.sqrt(sum([anchoring[i + 3] ** 2 for i in range(3)])) + oh = np.sqrt(sum([anchoring[i + 6] ** 2 for i in range(3)])) xspacing = int(width * grid_spacing / ow) yspacing = int(height * grid_spacing / oh) return xspacing, yspacing + def create_damage_mask(section, grid_spacing): width = section["width"] height = section["height"] @@ -65,7 +67,10 @@ def create_damage_mask(section, grid_spacing): y_coords = np.arange(gridy, height, yspacing) num_markers = len(grid_values) - markers = [(x_coords[i % len(x_coords)], y_coords[i // len(x_coords)]) for i in range(num_markers)] + markers = [ + (x_coords[i % len(x_coords)], y_coords[i // len(x_coords)]) + for i in range(num_markers) + ] binary_image = np.ones((len(y_coords), len(x_coords)), dtype=int) @@ -75,6 +80,7 @@ def create_damage_mask(section, grid_spacing): return binary_image + def folder_to_atlas_space( folder, quint_alignment, @@ -104,9 +110,15 @@ def folder_to_atlas_space( slices, gridspacing = load_visualign_json(quint_alignment) segmentations = get_segmentations(folder) flat_files, flat_file_nrs = get_flat_files(folder, use_flat) - points_list, centroids_list, region_areas_list, centroids_labels, points_labels,per_point_undamaged_list,per_centroid_undamaged_list = ( - initialize_lists(len(segmentations)) - ) + ( + points_list, + centroids_list, + region_areas_list, + centroids_labels, + points_labels, + per_point_undamaged_list, + per_centroid_undamaged_list, + ) = initialize_lists(len(segmentations)) threads = create_threads( segmentations, slices, @@ -125,11 +137,25 @@ def folder_to_atlas_space( object_cutoff, atlas_volume, use_flat, - gridspacing + gridspacing, ) start_and_join_threads(threads) - points, centroids, points_labels, centroids_labels, points_len, centroids_len, per_point_undamaged_list, per_centroid_undamaged_list = ( - process_results(points_list, centroids_list, points_labels, centroids_labels, per_point_undamaged_list, per_centroid_undamaged_list) + ( + points, + centroids, + points_labels, + centroids_labels, + points_len, + centroids_len, + per_point_undamaged_list, + per_centroid_undamaged_list, + ) = process_results( + points_list, + centroids_list, + points_labels, + centroids_labels, + per_point_undamaged_list, + per_centroid_undamaged_list, ) return ( points, @@ -141,7 +167,7 @@ def folder_to_atlas_space( centroids_len, segmentations, per_point_undamaged_list, - per_centroid_undamaged_list + per_centroid_undamaged_list, ) @@ -183,7 +209,7 @@ def initialize_lists(length): centroids_labels, points_labels, per_point_undamaged_list, - per_centroid_undamaged_list + per_centroid_undamaged_list, ) @@ -205,7 +231,7 @@ def create_threads( object_cutoff, atlas_volume, use_flat, - gridspacing + gridspacing, ): """ Creates threads for processing segmentations. @@ -261,7 +287,7 @@ def create_threads( object_cutoff, atlas_volume, use_flat, - gridspacing + gridspacing, ), ) threads.append(x) @@ -309,7 +335,7 @@ def get_region_areas( slice_dict, atlas_volume, triangulation, - damage_mask + damage_mask, ): """ Gets the region areas. @@ -327,11 +353,16 @@ def get_region_areas( Returns: DataFrame: DataFrame with region areas. """ - atlas_map = load_image(flat_file_atlas,slice_dict["anchoring"], atlas_volume, triangulation, (seg_width, seg_height), atlas_labels) - - region_areas = flat_to_dataframe( - atlas_map, damage_mask, (seg_width, seg_height) + atlas_map = load_image( + flat_file_atlas, + slice_dict["anchoring"], + atlas_volume, + triangulation, + (seg_width, seg_height), + atlas_labels, ) + + region_areas = flat_to_dataframe(atlas_map, damage_mask, (seg_width, seg_height)) return region_areas, atlas_map @@ -392,7 +423,7 @@ def segmentation_to_atlas_space( slice_dict, atlas_volume, triangulation, - damage_mask + damage_mask, ) atlas_map = rescale_image(atlas_map, (reg_width, reg_height)) y_scale, x_scale = transform_to_registration( @@ -411,16 +442,21 @@ def segmentation_to_atlas_space( np.round(scaled_centroidsX).astype(int), np.round(scaled_centroidsY).astype(int) ] if damage_mask is not None: - damage_mask = cv2.resize(damage_mask.astype(np.uint8), (atlas_map.shape[::-1]), interpolation=cv2.INTER_NEAREST).astype(bool) + damage_mask = cv2.resize( + damage_mask.astype(np.uint8), + (atlas_map.shape[::-1]), + interpolation=cv2.INTER_NEAREST, + ).astype(bool) per_point_undamaged = damage_mask[ - np.round(scaled_x).astype(int), np.round(scaled_y).astype(int) - ] + np.round(scaled_x).astype(int), np.round(scaled_y).astype(int) + ] per_centroid_undamaged = damage_mask[ - np.round(scaled_centroidsX).astype(int), np.round(scaled_centroidsY).astype(int) - ] + np.round(scaled_centroidsX).astype(int), + np.round(scaled_centroidsY).astype(int), + ] else: per_point_undamaged = np.ones(scaled_x.shape, dtype=bool) - per_centroid_undamaged = np.ones(scaled_centroidsX.shape, dtype=bool) + per_centroid_undamaged = np.ones(scaled_centroidsX.shape, dtype=bool) per_point_labels = per_point_labels[per_point_undamaged] per_centroid_labels = per_centroid_labels[per_centroid_undamaged] new_x, new_y, centroids_new_x, centroids_new_y = get_transformed_coordinates( @@ -451,7 +487,9 @@ def segmentation_to_atlas_space( per_centroid_undamaged if centroids is not None else [] ) points_labels[index] = np.array(per_point_labels if points is not None else []) - per_point_undamaged_list[index] = np.array(per_point_undamaged if points is not None else []) + per_point_undamaged_list[index] = np.array( + per_point_undamaged if points is not None else [] + ) def get_triangulation(slice_dict, reg_width, reg_height, non_linear): diff --git a/PyNutil/processing/counting_and_load.py b/PyNutil/processing/counting_and_load.py index da1fce6f65091054f932022d4620b6c5ba45fe24..0b549e6a4037cdb5ebe0e0e021d969376c47b8a1 100644 --- a/PyNutil/processing/counting_and_load.py +++ b/PyNutil/processing/counting_and_load.py @@ -5,6 +5,7 @@ import cv2 from .generate_target_slice import generate_target_slice from .visualign_deformations import transform_vec + # related to counting and load def label_points(points, label_volume, scale_factor=1): """ @@ -49,7 +50,11 @@ def label_points(points, label_volume, scale_factor=1): # related to counting_and_load def pixel_count_per_region( - labels_dict_points, labeled_dict_centroids, current_points_undamaged, current_centroids_undamaged, df_label_colours + labels_dict_points, + labeled_dict_centroids, + current_points_undamaged, + current_centroids_undamaged, + df_label_colours, ): """ Counts the number of pixels per region and writes to a DataFrame. @@ -77,75 +82,59 @@ def pixel_count_per_region( # Which regions have pixels, and how many pixels are there per region # Create a list of unique regions and pixel counts per region counts_per_label = { - "idx" : [], - "name" : [], - "r" : [], - "g" : [], - "b" : [], - "pixel_count" : [], - "undamaged_pixel_count" : [], - "damaged_pixel_counts" : [], - "object_count" : [], - "undamaged_object_count" : [], - "damaged_object_count" : [], - } + "idx": [], + "name": [], + "r": [], + "g": [], + "b": [], + "pixel_count": [], + "undamaged_pixel_count": [], + "damaged_pixel_counts": [], + "object_count": [], + "undamaged_object_count": [], + "damaged_object_count": [], + } for index, row in df_label_colours.iterrows(): if row["idx"] in counted_labels_points_undamaged: - clpu = label_counts_points_undamaged[counted_labels_points_undamaged == row["idx"]][0] + clpu = label_counts_points_undamaged[ + counted_labels_points_undamaged == row["idx"] + ][0] else: clpu = 0 if row["idx"] in counted_labels_points_damaged: - clpd = label_counts_points_damaged[counted_labels_points_damaged == row["idx"]][0] + clpd = label_counts_points_damaged[ + counted_labels_points_damaged == row["idx"] + ][0] else: clpd = 0 if row["idx"] in counted_labels_centroids_undamaged: - clcu = counted_labels_centroids_undamaged[counted_labels_centroids_undamaged == row["idx"]][0] + clcu = counted_labels_centroids_undamaged[ + counted_labels_centroids_undamaged == row["idx"] + ][0] else: clcu = 0 if row["idx"] in counted_labels_centroids_damaged: - clcd = counted_labels_centroids_damaged[counted_labels_centroids_damaged == row["idx"]][0] + clcd = counted_labels_centroids_damaged[ + counted_labels_centroids_damaged == row["idx"] + ][0] else: clcd = 0 - if clcd==clcu==clpd==clpu==0: + if clcd == clcu == clpd == clpu == 0: continue - counts_per_label["idx"].append( - row["idx"] - ) - counts_per_label["name"].append( - row["name"] - ) - counts_per_label["r"].append( - int(row["r"]) - ) - counts_per_label["g"].append( - int(row["g"]) - ) - counts_per_label["b"].append( - int(row["b"]) - ) - counts_per_label["pixel_count"].append( - clpu + clpd - ) - counts_per_label["undamaged_pixel_count"].append( - clpu - ) - counts_per_label["damaged_pixel_counts"].append( - clpd - ) - counts_per_label["object_count"].append( - clcu + clcd - ) - counts_per_label["undamaged_object_count"].append( - clcu - ) - counts_per_label["damaged_object_count"].append( - clcd - ) - - df_counts_per_label = pd.DataFrame( - counts_per_label - ) + counts_per_label["idx"].append(row["idx"]) + counts_per_label["name"].append(row["name"]) + counts_per_label["r"].append(int(row["r"])) + counts_per_label["g"].append(int(row["g"])) + counts_per_label["b"].append(int(row["b"])) + counts_per_label["pixel_count"].append(clpu + clpd) + counts_per_label["undamaged_pixel_count"].append(clpu) + counts_per_label["damaged_pixel_counts"].append(clpd) + counts_per_label["object_count"].append(clcu + clcd) + counts_per_label["undamaged_object_count"].append(clcu) + counts_per_label["damaged_object_count"].append(clcd) + + df_counts_per_label = pd.DataFrame(counts_per_label) return df_counts_per_label @@ -311,10 +300,7 @@ def warp_image(image, triangulation, rescaleXY): return new_image -def flat_to_dataframe( - image, - damage_mask, - rescaleXY=None): +def flat_to_dataframe(image, damage_mask, rescaleXY=None): """ Converts a flat file to a DataFrame. @@ -332,13 +318,33 @@ def flat_to_dataframe( """ scale_factor = calculate_scale_factor(image, rescaleXY) if damage_mask is not None: - damage_mask = cv2.resize(damage_mask.astype(np.uint8), (image.shape[::-1]), interpolation=cv2.INTER_NEAREST).astype(bool) - undamaged_df_area_per_label = count_pixels_per_label(image[damage_mask], scale_factor) - damaged_df_area_per_label = count_pixels_per_label(image[~damage_mask], scale_factor) - undamaged_df_area_per_label = undamaged_df_area_per_label.rename(columns={"region_area": "undamaged_region_area"}) - damaged_df_area_per_label = damaged_df_area_per_label.rename(columns={"region_area": "damaged_region_area"}) - df_area_per_label = pd.merge(undamaged_df_area_per_label, damaged_df_area_per_label, on='idx', how='outer').fillna(0) - df_area_per_label["region_area"] = df_area_per_label["undamaged_region_area"] + df_area_per_label["damaged_region_area"] + damage_mask = cv2.resize( + damage_mask.astype(np.uint8), + (image.shape[::-1]), + interpolation=cv2.INTER_NEAREST, + ).astype(bool) + undamaged_df_area_per_label = count_pixels_per_label( + image[damage_mask], scale_factor + ) + damaged_df_area_per_label = count_pixels_per_label( + image[~damage_mask], scale_factor + ) + undamaged_df_area_per_label = undamaged_df_area_per_label.rename( + columns={"region_area": "undamaged_region_area"} + ) + damaged_df_area_per_label = damaged_df_area_per_label.rename( + columns={"region_area": "damaged_region_area"} + ) + df_area_per_label = pd.merge( + undamaged_df_area_per_label, + damaged_df_area_per_label, + on="idx", + how="outer", + ).fillna(0) + df_area_per_label["region_area"] = ( + df_area_per_label["undamaged_region_area"] + + df_area_per_label["damaged_region_area"] + ) else: df_area_per_label = count_pixels_per_label(image, scale_factor) df_area_per_label["undamaged_region_area"] = df_area_per_label["region_area"] diff --git a/PyNutil/processing/data_analysis.py b/PyNutil/processing/data_analysis.py index 64cadce6b9942512fdd0a3366b7e61719637532e..dc7816afdebc3dcc2c3972145ee9025206d8ef30 100644 --- a/PyNutil/processing/data_analysis.py +++ b/PyNutil/processing/data_analysis.py @@ -2,6 +2,7 @@ import pandas as pd from .counting_and_load import pixel_count_per_region import numpy as np + def map_to_custom_regions(custom_regions_dict, points_labels): custom_points_labels = np.zeros_like(points_labels) for i in np.unique(points_labels): @@ -74,7 +75,9 @@ def apply_custom_regions(df, custom_regions_dict): grouped_df = grouped_df.rename(columns={"custom_region_name": "name"}) grouped_df["area_fraction"] = grouped_df["pixel_count"] / grouped_df["region_area"] - grouped_df["undamaged_area_fraction"] = grouped_df["undamaged_pixel_count"] / grouped_df["undamaged_region_area"] + grouped_df["undamaged_area_fraction"] = ( + grouped_df["undamaged_pixel_count"] / grouped_df["undamaged_region_area"] + ) common_columns = [col for col in df.columns if col in grouped_df.columns] grouped_df = grouped_df.reindex( columns=common_columns @@ -82,6 +85,7 @@ def apply_custom_regions(df, custom_regions_dict): ) return grouped_df, df + def quantify_labeled_points( points_len, centroids_len, @@ -90,7 +94,7 @@ def quantify_labeled_points( labeled_points_centroids, atlas_labels, per_point_undamaged, - per_centroid_undamaged + per_centroid_undamaged, ): """ Quantifies labeled points and returns various DataFrames. @@ -118,7 +122,7 @@ def quantify_labeled_points( region_areas_list, atlas_labels, per_point_undamaged, - per_centroid_undamaged + per_centroid_undamaged, ) label_df = _combine_slice_reports(per_section_df, atlas_labels) @@ -133,7 +137,7 @@ def _quantify_per_section( region_areas_list, atlas_labels, per_point_undamaged, - per_centroid_undamaged + per_centroid_undamaged, ): """ Quantifies labeled points per section. @@ -159,7 +163,11 @@ def _quantify_per_section( current_points_undamaged = per_point_undamaged[prev_pl : prev_pl + pl] current_centroids_undamaged = per_centroid_undamaged[prev_cl : prev_cl + cl] current_df = pixel_count_per_region( - current_points, current_centroids, current_points_undamaged, current_centroids_undamaged, atlas_labels + current_points, + current_centroids, + current_points_undamaged, + current_centroids_undamaged, + atlas_labels, ) current_df_new = _merge_dataframes(current_df, ra, atlas_labels) per_section_df.append(current_df_new) @@ -211,7 +219,9 @@ def _combine_slice_reports(per_section_df, atlas_labels): .drop(columns=["area_fraction"]) ) label_df["area_fraction"] = label_df["pixel_count"] / label_df["region_area"] - label_df["undamaged_area_fraction"] = label_df["undamaged_pixel_count"] / label_df["undamaged_region_area"] + label_df["undamaged_area_fraction"] = ( + label_df["undamaged_pixel_count"] / label_df["undamaged_region_area"] + ) label_df.fillna(0, inplace=True) label_df = label_df.set_index("idx") diff --git a/PyNutil/processing/generate_target_slice.py b/PyNutil/processing/generate_target_slice.py index 937d4d474047f7c8484744c2e50c9d11a47c58b1..a8bc9c570d4d43b7d4db49e0f73c1cfd44bfb686 100644 --- a/PyNutil/processing/generate_target_slice.py +++ b/PyNutil/processing/generate_target_slice.py @@ -1,6 +1,7 @@ import numpy as np import math + def generate_target_slice(ouv, atlas): """ Generates a target slice from the atlas using the given orientation and vectors. diff --git a/PyNutil/processing/transformations.py b/PyNutil/processing/transformations.py index 671ea2d5e9789bac7a6d2e5210a480e0a3096a1d..62b8ddb4ebdc8b18ce4fb05c07b45fd65114f0e4 100644 --- a/PyNutil/processing/transformations.py +++ b/PyNutil/processing/transformations.py @@ -1,6 +1,7 @@ import numpy as np from .visualign_deformations import transform_vec + def transform_to_registration(seg_height, seg_width, reg_height, reg_width): """ Returns the scaling factors to transform the segmentation to the registration space. diff --git a/PyNutil/processing/utils.py b/PyNutil/processing/utils.py index b09bc8323834a577b6513e04af1a24e34b82fcbe..ed9b64c554e733fd9675420a69cfc66aa46ced61 100644 --- a/PyNutil/processing/utils.py +++ b/PyNutil/processing/utils.py @@ -3,6 +3,7 @@ import pandas as pd import re from glob import glob + def number_sections(filenames, legacy=False): """ Returns the section numbers of filenames. @@ -196,7 +197,14 @@ def start_and_join_threads(threads): [t.join() for t in threads] -def process_results(points_list, centroids_list, points_labels, centroids_labels, points_undamaged_list, centroids_undamaged_list): +def process_results( + points_list, + centroids_list, + points_labels, + centroids_labels, + points_undamaged_list, + centroids_undamaged_list, +): """ Processes the results from the threads. @@ -218,7 +226,9 @@ def process_results(points_list, centroids_list, points_labels, centroids_labels points_labels = [pl for pl in points_labels if None not in pl] centroids_labels = [cl for cl in centroids_labels if None not in cl] points_undamaged_list = [pul for pul in points_undamaged_list if None not in pul] - centroids_undamaged_list = [cul for cul in centroids_undamaged_list if None not in cul] + centroids_undamaged_list = [ + cul for cul in centroids_undamaged_list if None not in cul + ] if len(points_list) == 0: points = np.array([]) points_labels = np.array([]) @@ -237,4 +247,13 @@ def process_results(points_list, centroids_list, points_labels, centroids_labels centroids_labels = np.concatenate(centroids_labels) centroids_undamaged = np.concatenate(centroids_undamaged_list) - return points, centroids, points_labels, centroids_labels, points_len, centroids_len, points_undamaged, centroids_undamaged + return ( + points, + centroids, + points_labels, + centroids_labels, + points_len, + centroids_len, + points_undamaged, + centroids_undamaged, + ) diff --git a/PyNutil/processing/visualign_deformations.py b/PyNutil/processing/visualign_deformations.py index 3b6e12ff5c88736aca1801d1dd167c6fef1c069c..d0303452215f11f5fc399124621e02cf6cb91759 100644 --- a/PyNutil/processing/visualign_deformations.py +++ b/PyNutil/processing/visualign_deformations.py @@ -1,6 +1,8 @@ """This code was written by Gergely Csucs, Harry Carey and Rembrandt Bakker""" + import numpy as np + def triangulate(w, h, markers): """ Triangulates a set of markers. diff --git a/demos/basic_example_custom_atlas.py b/demos/basic_example_custom_atlas.py index b2b9c189f5e61108bd56a1c9037d2a3cc13ea78e..d0f4978d1beb77a768ec7557abd5ddf0390a02e2 100644 --- a/demos/basic_example_custom_atlas.py +++ b/demos/basic_example_custom_atlas.py @@ -1,4 +1,5 @@ import sys + sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from PyNutil import PyNutil import os