Add region areas to reports

Working session with Harry.
Read in flat files from a folder. Scale flat files to size of seg. Calculate region areas per slice. Create slice reports with region areas and object areas. Create report for all slices by summing slice reports. In progress.

PyNutil/coordinate_extraction.py

@@ -3,7 +3,7 @@ import pandas as pd
 from DeepSlice.coord_post_processing.spacing_and_indexing import number_sections
 import json
 from .read_and_write import load_visualign_json
-from .counting_and_load import label_points
+from .counting_and_load import label_points, flat_to_dataframe
 from .visualign_deformations import triangulate, transform_vec
 from glob import glob
 from tqdm import tqdm
@@ -87,6 +87,7 @@ def transform_to_atlas_space(anchoring, y, x, reg_height, reg_width):
 def folder_to_atlas_space(
     folder,
     quint_alignment,
+    atlas_labels,
     pixel_id=[0, 0, 0],
     non_linear=True,
     method="all",
@@ -103,25 +104,37 @@ def folder_to_atlas_space(
         for file in glob(folder + "/*")
         if any([file.endswith(type) for type in segmentation_file_types])
     ]
+    flat_files = [
+        file
+        for file in glob(folder + "/flat_files/*")
+        if any([file.endswith('.flat')])
+    ]
     # Order segmentations and section_numbers
     # segmentations = [x for _,x in sorted(zip(section_numbers,segmentations))]
     # section_numbers.sort()
     points_list = [None] * len(segmentations)
     centroids_list = [None] * len(segmentations)
+    region_areas_list = [None] * len(segmentations)
     threads = []
+    flat_file_nrs = [int(number_sections([ff])[0]) for ff in flat_files]
     for segmentation_path, index in zip(segmentations, range(len(segmentations))):
         seg_nr = int(number_sections([segmentation_path])[0])
         current_slice_index = np.where([s["nr"] == seg_nr for s in slices])
         current_slice = slices[current_slice_index[0][0]]
+        current_flat_file_index = np.where([f == seg_nr for f in flat_file_nrs])
+        current_flat = flat_files[current_flat_file_index[0][0]]
         x = threading.Thread(
             target=segmentation_to_atlas_space,
             args=(
                 current_slice,
                 segmentation_path,
+                atlas_labels, 
+                current_flat,
                 pixel_id,
                 non_linear,
                 points_list,
                 centroids_list,
+                region_areas_list,
                 index,
                 method,
                 object_cutoff,
@@ -143,6 +156,7 @@ def folder_to_atlas_space(
     return (
         np.array(points),
         np.array(centroids),
+        region_areas_list,
         points_len,
         centroids_len,
         segmentations,
@@ -152,13 +166,17 @@ def folder_to_atlas_space(
 def segmentation_to_atlas_space(
     slice,
     segmentation_path,
+    atlas_labels, 
+    flat_file_atlas,
     pixel_id="auto",
     non_linear=True,
     points_list=None,
     centroids_list=None,
+    region_areas_list=None, 
     index=None,
     method="per_pixel",
     object_cutoff=0,
+    
 ):
     """Combines many functions to convert a segmentation to atlas space. It takes care
     of deformations."""
@@ -177,6 +195,7 @@ def segmentation_to_atlas_space(
     seg_width = segmentation.shape[1]
     reg_height = slice["height"]
     reg_width = slice["width"]
+    region_areas = flat_to_dataframe(flat_file_atlas, atlas_labels, (seg_width,seg_height))
     # This calculates reg/seg
     y_scale, x_scale = transform_to_registration(
         seg_height, seg_width, reg_height, reg_width
@@ -226,6 +245,7 @@ def segmentation_to_atlas_space(
     )
     points_list[index] = np.array(points)
     centroids_list[index] = np.array(centroids)
+    region_areas_list[index] = region_areas
 
 
 def get_centroids(segmentation, pixel_id, y_scale, x_scale, object_cutoff=0):

PyNutil/counting_and_load.py

@@ -4,6 +4,7 @@ import struct
 import matplotlib.pyplot as plt
 import os
 import nrrd
+import cv2
 
 
 # related to counting and load
@@ -87,6 +88,7 @@ def pixel_count_per_region(
     # Find colours corresponding to each region ID and add to the pandas dataframe
 
     # Look up name, r, g, b in df_allen_colours in df_counts_per_label based on "idx"
+    #Sharon, remove this here
     new_rows = []
     for index, row in df_counts_per_label.iterrows():
         mask = df_label_colours["idx"] == row["idx"]
@@ -117,11 +119,14 @@ def pixel_count_per_region(
 """Read flat file and write into an np array"""
 """Read flat file, write into an np array, assign label file values, return array"""
 
-def flat_to_dataframe(flat_file, labelfile):
+def flat_to_dataframe(flat_file, labelfile, rescaleXY=False):
+    if rescaleXY:
+        rescaleX, rescaleY = rescaleXY
     with open(flat_file, "rb") as f:
         # I don't know what b is, w and h are the width and height that we get from the
         # flat file header
         b, w, h = struct.unpack(">BII", f.read(9))
+        print(w, h)
         # Data is a one dimensional list of values
         # It has the shape width times height
         data = struct.unpack(">" + ("xBH"[b] * (w * h)), f.read(b * w * h))
@@ -137,12 +142,15 @@ def flat_to_dataframe(flat_file, labelfile):
 
     image_arr = np.array(image)
     #return image_arr
+    if rescaleXY:
+        w,h= rescaleXY
+        image_arr = cv2.resize(image_arr, (h, w), interpolation=cv2.INTER_NEAREST)
+
 
     """assign label file values into image array"""
-    labelfile = pd.read_csv(labelfile)
     allen_id_image = np.zeros((h, w))  # create an empty image array
     coordsy, coordsx = np.meshgrid(list(range(w)), list(range(h)))
-    values = image_arr[coordsx, coordsy]  # assign x,y coords from image_array into values
+    values = image_arr[coordsy,coordsx]  # assign x,y coords from image_array into values
     lbidx = labelfile["idx"].values
     allen_id_image = lbidx[values.astype(int)]
     #return allen_id_image

PyNutil/load_workflow.py

@@ -3,6 +3,7 @@
 import struct
 import matplotlib.pyplot as plt
 import pandas as pd
+import cv2
 import numpy as np
 
 #from read_and_write import flat_to_array, label_to_array
@@ -10,12 +11,16 @@ from counting_and_load import flat_to_dataframe
 
 base = r"../test_data/tTA_2877_NOP_s037_atlasmap/2877_NOP_tTA_lacZ_Xgal_s037_nl.flat"
 label = r"../annotation_volumes\allen2017_colours.csv"
-
+##optional
+seg = r"../test_data/tTA_2877_NOP_s037_seg/2877_NOP_tTA_lacZ_Xgal_resize_Simple_Seg_s037.png"
+segim = cv2.imread(seg)
+#the indexing [:2] means the first two values and [::-1] means reverse the list
+segXY = segim.shape[:2][::-1]
 #image_arr = flat_to_array(base, label)
 
 #plt.imshow(flat_to_array(base, label))
 
-df_area_per_label = flat_to_dataframe(base, label)
+df_area_per_label = flat_to_dataframe(base, label, segXY)
 
 """count pixels in np array for unique idx, return pd df"""
 #unique_ids, counts = np.unique(allen_id_image, return_counts=True)
@@ -54,7 +59,7 @@ df_area_per_label_name = pd.DataFrame(new_rows)
 
 print(df_area_per_label_name)
 df_area_per_label_name.to_csv(
-    "../outputs/test5_s037_area_per_idx.csv", sep=";", na_rep="", index=False
+    "../outputs/NOP_s037_regionareas.csv", sep=";", na_rep="", index=False
 )
 
 

PyNutil/main.py

@@ -163,12 +163,14 @@ class PyNutil:
         (
             pixel_points,
             centroids,
+            region_areas_list,
             points_len,
             centroids_len,
             segmentation_filenames,
         ) = folder_to_atlas_space(
             self.segmentation_folder,
             self.alignment_json,
+            self.atlas_labels,
             pixel_id=self.colour,
             non_linear=non_linear,
             method=method,
@@ -181,6 +183,7 @@ class PyNutil:
         self.points_len = points_len
         self.centroids_len = centroids_len
         self.segmentation_filenames = segmentation_filenames
+        self.region_areas_list = region_areas_list
 
     def quantify_coordinates(self):
         """Quantifies the pixel coordinates by region.
@@ -207,15 +210,12 @@ class PyNutil:
                 self.pixel_points, self.atlas_volume, scale_factor=1
             )
 
-        self.label_df = pixel_count_per_region(
-            labeled_points, labeled_points_centroids, self.atlas_labels
-        )
         prev_pl = 0
         prev_cl = 0
         per_section_df = []
         current_centroids = None
         current_points = None
-        for pl,cl in zip(self.points_len, self.centroids_len):
+        for pl,cl,ra in zip(self.points_len, self.centroids_len, self.region_areas_list):
             if hasattr(self, "centroids"):
                 current_centroids = labeled_points_centroids[prev_cl : prev_cl + cl]
             if hasattr(self, "pixel_points"):
@@ -223,9 +223,14 @@ class PyNutil:
             current_df = pixel_count_per_region(
                 current_points, current_centroids, self.atlas_labels
             )
+
+            current_df = ra.merge(current_df, on='idx', how='left')
+            ##Sharon. I would guess you should add the rgb and name adding code here
             per_section_df.append(current_df)
             prev_pl += pl
             prev_cl += cl
+        ##Sharon. and then here you should group on r,g,b,idx, and name since you dont want any of these summed
+        self.label_df =  pd.concat(per_section_df).groupby(['idx']).sum().reset_index()
 
         self.labeled_points = labeled_points
         self.labeled_points_centroids = labeled_points_centroids