From b1517a4755f9b9bf2c6209f0b0df4739f4a8b61a Mon Sep 17 00:00:00 2001
From: Harry Carey <harry.carey95@gmail.com>
Date: Fri, 13 Oct 2023 15:48:39 +0200
Subject: [PATCH] black reformatting
---
PyNutil/coordinate_extraction.py | 75 +++---
PyNutil/counting_and_load.py | 59 +++--
PyNutil/generate_target_slice.py | 8 +-
PyNutil/load_workflow.py | 14 +-
PyNutil/main.py | 67 +++---
PyNutil/propagation.py | 231 ++++++++++---------
PyNutil/read_and_write.py | 50 ++--
PyNutil/reconstruct_dzi.py | 15 +-
PyNutil/visualign_deformations.py | 14 +-
messing_around_files/pyflat.py | 40 ++--
messing_around_files/testing_openflatfile.py | 10 +-
scripts/create_test_data.py | 29 ++-
testOOP.py | 4 +-
13 files changed, 353 insertions(+), 263 deletions(-)
diff --git a/PyNutil/coordinate_extraction.py b/PyNutil/coordinate_extraction.py
index 6079f2c..157bdc6 100644
--- a/PyNutil/coordinate_extraction.py
+++ b/PyNutil/coordinate_extraction.py
@@ -9,7 +9,7 @@ from tqdm import tqdm
import cv2
from skimage import measure
import threading
-import re
+import re
from .reconstruct_dzi import reconstruct_dzi
@@ -42,6 +42,7 @@ def number_sections(filenames, legacy=False):
raise ValueError("No section numbers found in filenames")
return section_numbers
+
# related to coordinate_extraction
def get_centroids_and_area(segmentation, pixel_cut_off=0):
"""This function returns the center coordinate of each object in the segmentation.
@@ -122,7 +123,7 @@ def folder_to_atlas_space(
method="all",
object_cutoff=0,
atlas_volume=None,
- use_flat = False,
+ use_flat=False,
):
"""Apply Segmentation to atlas space to all segmentations in a folder."""
"""Return pixel_points, centroids, points_len, centroids_len, segmentation_filenames, """
@@ -144,29 +145,30 @@ def folder_to_atlas_space(
flat_files = [
file
for file in glob(folder + "/flat_files/*")
- if any([file.endswith('.flat'), file.endswith('.seg')])
+ if any([file.endswith(".flat"), file.endswith(".seg")])
]
print(f"Found {len(flat_files)} flat files in folder {folder}")
flat_file_nrs = [int(number_sections([ff])[0]) for ff in flat_files]
-
# Order segmentations and section_numbers
# segmentations = [x for _,x in sorted(zip(section_numbers,segmentations))]
# section_numbers.sort()
points_list = [np.array([])] * len(segmentations)
centroids_list = [np.array([])] * len(segmentations)
region_areas_list = [
- pd.DataFrame({
- "idx": [],
- "name": [],
- "r": [],
- "g": [],
- "b": [],
- "region_area": [],
- "pixel_count": [],
- "object_count": [],
- "area_fraction": []
- })
+ pd.DataFrame(
+ {
+ "idx": [],
+ "name": [],
+ "r": [],
+ "g": [],
+ "b": [],
+ "region_area": [],
+ "pixel_count": [],
+ "object_count": [],
+ "area_fraction": [],
+ }
+ )
] * len(segmentations)
threads = []
for segmentation_path, index in zip(segmentations, range(len(segmentations))):
@@ -186,7 +188,7 @@ def folder_to_atlas_space(
args=(
current_slice,
segmentation_path,
- atlas_labels,
+ atlas_labels,
current_flat,
pixel_id,
non_linear,
@@ -197,7 +199,7 @@ def folder_to_atlas_space(
method,
object_cutoff,
atlas_volume,
- use_flat
+ use_flat,
),
)
threads.append(x)
@@ -207,14 +209,14 @@ def folder_to_atlas_space(
# Wait for threads to finish
[t.join() for t in threads]
# Flatten points_list
-
points_len = [len(points) if points is not None else 0 for points in points_list]
- centroids_len = [len(centroids) if centroids is not None else 0 for centroids in centroids_list]
+ centroids_len = [
+ len(centroids) if centroids is not None else 0 for centroids in centroids_list
+ ]
points = np.concatenate(points_list)
centroids = np.concatenate(centroids_list)
-
return (
np.array(points),
np.array(centroids),
@@ -228,19 +230,18 @@ def folder_to_atlas_space(
def segmentation_to_atlas_space(
slice,
segmentation_path,
- atlas_labels,
+ atlas_labels,
flat_file_atlas=None,
pixel_id="auto",
non_linear=True,
points_list=None,
centroids_list=None,
- region_areas_list=None,
+ region_areas_list=None,
index=None,
method="per_pixel",
object_cutoff=0,
atlas_volume=None,
- use_flat=False
-
+ use_flat=False,
):
"""Combines many functions to convert a segmentation to atlas space. It takes care
of deformations."""
@@ -261,7 +262,7 @@ def segmentation_to_atlas_space(
# Currently only works for a single label
print("length of non background pixels: ", len(segmentation_no_background))
pixel_id = segmentation_no_background[0]
- print('detected pixel_id: ', pixel_id)
+ print("detected pixel_id: ", pixel_id)
# Transform pixels to registration space (the registered image and segmentation have different dimensions)
seg_height = segmentation.shape[0]
@@ -269,9 +270,17 @@ def segmentation_to_atlas_space(
reg_height = slice["height"]
reg_width = slice["width"]
if use_flat == True:
- region_areas = flat_to_dataframe(atlas_labels, flat_file_atlas, (seg_width,seg_height))
+ region_areas = flat_to_dataframe(
+ atlas_labels, flat_file_atlas, (seg_width, seg_height)
+ )
else:
- region_areas = flat_to_dataframe(atlas_labels, flat_file_atlas, (seg_width,seg_height), slice["anchoring"], atlas_volume)
+ region_areas = flat_to_dataframe(
+ atlas_labels,
+ flat_file_atlas,
+ (seg_width, seg_height),
+ slice["anchoring"],
+ atlas_volume,
+ )
# This calculates reg/seg
y_scale, x_scale = transform_to_registration(
seg_height, seg_width, reg_height, reg_width
@@ -287,7 +296,6 @@ def segmentation_to_atlas_space(
if non_linear:
if "markers" in slice:
-
# This creates a triangulation using the reg width
triangulation = triangulate(reg_width, reg_height, slice["markers"])
if method in ["per_pixel", "all"]:
@@ -301,7 +309,10 @@ def segmentation_to_atlas_space(
triangulation, scaled_centroidsX, scaled_centroidsY
)
else:
- centroids_new_x, centroids_new_y = scaled_centroidsX, scaled_centroidsY
+ centroids_new_x, centroids_new_y = (
+ scaled_centroidsX,
+ scaled_centroidsY,
+ )
else:
print(
f"No markers found for {slice['filename']}, result for section will be linear."
@@ -326,7 +337,11 @@ def segmentation_to_atlas_space(
if method in ["per_object", "all"]:
if centroids_new_x is not None:
centroids = transform_to_atlas_space(
- slice["anchoring"], centroids_new_y, centroids_new_x, reg_height, reg_width
+ slice["anchoring"],
+ centroids_new_y,
+ centroids_new_x,
+ reg_height,
+ reg_width,
)
else:
centroids = np.array([])
diff --git a/PyNutil/counting_and_load.py b/PyNutil/counting_and_load.py
index 9dce6ec..9da941f 100644
--- a/PyNutil/counting_and_load.py
+++ b/PyNutil/counting_and_load.py
@@ -89,7 +89,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
+ # Sharon, remove this here
new_rows = []
for index, row in df_counts_per_label.iterrows():
mask = df_label_colours["idx"] == row["idx"]
@@ -106,12 +106,14 @@ def pixel_count_per_region(
new_rows.append(row)
- df_counts_per_label_name = pd.DataFrame(new_rows, columns=["idx", "name", "pixel_count", "object_count", "r", "g", "b"])
- #Task for Sharon:
- #If you can get the areas per region from the flat file here
- #you can then use those areas to calculate the load per region here
+ df_counts_per_label_name = pd.DataFrame(
+ new_rows, columns=["idx", "name", "pixel_count", "object_count", "r", "g", "b"]
+ )
+ # Task for Sharon:
+ # If you can get the areas per region from the flat file here
+ # you can then use those areas to calculate the load per region here
# and add to dataframe
- #see messing around pyflat.py
+ # see messing around pyflat.py
return df_counts_per_label_name
@@ -123,6 +125,7 @@ import cv2
import numpy as np
import pandas as pd
+
def read_flat_file(file):
with open(file, "rb") as f:
b, w, h = struct.unpack(">BII", f.read(9))
@@ -134,51 +137,61 @@ def read_flat_file(file):
image[y, x] = image_data[x + y * w]
return image
+
def read_seg_file(file):
- with open(file,"rb") as f:
+ with open(file, "rb") as f:
+
def byte():
return f.read(1)[0]
+
def code():
- c=byte()
- if(c<0):
+ c = byte()
+ if c < 0:
raise "!"
- return c if c<128 else (c&127)|(code()<<7)
+ return c if c < 128 else (c & 127) | (code() << 7)
+
if "SegRLEv1" != f.read(8).decode():
raise "Header mismatch"
- atlas=f.read(code()).decode()
+ atlas = f.read(code()).decode()
print(f"Target atlas: {atlas}")
- codes=[code() for x in range(code())]
- w=code()
- h=code()
- data=[]
- while len(data)<w*h:
- data += [codes[byte() if len(codes)<=256 else code()]]*(code()+1)
+ codes = [code() for x in range(code())]
+ w = code()
+ h = code()
+ data = []
+ while len(data) < w * h:
+ data += [codes[byte() if len(codes) <= 256 else code()]] * (code() + 1)
image_data = np.array(data)
image = np.reshape(image_data, (h, w))
return image
+
def rescale_image(image, rescaleXY):
w, h = rescaleXY
return cv2.resize(image, (h, w), interpolation=cv2.INTER_NEAREST)
+
def assign_labels_to_image(image, labelfile):
- w,h = image.shape
+ w, h = image.shape
allen_id_image = np.zeros((h, w)) # create an empty image array
coordsy, coordsx = np.meshgrid(list(range(w)), list(range(h)))
-
+
values = image[coordsy, coordsx]
lbidx = labelfile["idx"].values
allen_id_image = lbidx[values.astype(int)]
return allen_id_image
+
def count_pixels_per_label(image):
unique_ids, counts = np.unique(image, return_counts=True)
area_per_label = list(zip(unique_ids, counts))
df_area_per_label = pd.DataFrame(area_per_label, columns=["idx", "region_area"])
return df_area_per_label
-def flat_to_dataframe(labelfile, file=None, rescaleXY=None, image_vector=None, volume=None):
+
+def flat_to_dataframe(
+ labelfile, file=None, rescaleXY=None, image_vector=None, volume=None
+):
if (image_vector is not None) and (volume is not None):
image = generate_target_slice(image_vector, volume)
image = np.float64(image)
@@ -186,8 +199,8 @@ def flat_to_dataframe(labelfile, file=None, rescaleXY=None, image_vector=None,
image = read_flat_file(file)
elif file.endswith(".seg"):
image = read_seg_file(file)
- print('datatype', image.dtype)
- print('image shape open', image.shape)
+ print("datatype", image.dtype)
+ print("image shape open", image.shape)
if rescaleXY:
image = rescale_image(image, rescaleXY)
@@ -196,4 +209,4 @@ def flat_to_dataframe(labelfile, file=None, rescaleXY=None, image_vector=None,
else:
allen_id_image = image
df_area_per_label = count_pixels_per_label(allen_id_image)
- return df_area_per_label
\ No newline at end of file
+ return df_area_per_label
diff --git a/PyNutil/generate_target_slice.py b/PyNutil/generate_target_slice.py
index 9cd83c7..31681b8 100644
--- a/PyNutil/generate_target_slice.py
+++ b/PyNutil/generate_target_slice.py
@@ -29,8 +29,12 @@ def generate_target_slice(alignment, volume):
Z_Coords = np.round(Z_Coords).astype(int)
out_bounds_Coords = (
- (X_Coords > bounds[0]) | (Y_Coords > bounds[1]) | (Z_Coords > bounds[2])
- | (X_Coords < 0) | (Y_Coords < 0) | (Z_Coords < 0)
+ (X_Coords > bounds[0])
+ | (Y_Coords > bounds[1])
+ | (Z_Coords > bounds[2])
+ | (X_Coords < 0)
+ | (Y_Coords < 0)
+ | (Z_Coords < 0)
)
X_pad = X_Coords.copy()
Y_pad = Y_Coords.copy()
diff --git a/PyNutil/load_workflow.py b/PyNutil/load_workflow.py
index ac3dd3d..9675196 100644
--- a/PyNutil/load_workflow.py
+++ b/PyNutil/load_workflow.py
@@ -6,7 +6,7 @@ import pandas as pd
import cv2
import numpy as np
-#from read_and_write import flat_to_array, label_to_array
+# from read_and_write import flat_to_array, label_to_array
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"
@@ -14,21 +14,21 @@ 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
+# 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)
+# image_arr = flat_to_array(base, label)
-#plt.imshow(flat_to_array(base, label))
+# plt.imshow(flat_to_array(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)
+# unique_ids, counts = np.unique(allen_id_image, return_counts=True)
-#area_per_label = list(zip(unique_ids, counts))
+# area_per_label = list(zip(unique_ids, counts))
# create a list of unique regions and pixel counts per region
-#df_area_per_label = pd.DataFrame(area_per_label, columns=["idx", "area_count"])
+# df_area_per_label = pd.DataFrame(area_per_label, columns=["idx", "area_count"])
# create a pandas df with regions and pixel counts
diff --git a/PyNutil/main.py b/PyNutil/main.py
index 6277a1f..44ac2d4 100644
--- a/PyNutil/main.py
+++ b/PyNutil/main.py
@@ -8,7 +8,6 @@ from datetime import datetime
import os
-
class PyNutil:
"""A utility class for working with brain atlases and segmentation data.
@@ -132,9 +131,10 @@ class PyNutil:
atlas_labels = pd.read_csv(f"{atlas_root_path}{atlas_label_path}")
print("atlas labels loaded ✅")
return atlas_volume, atlas_labels
-
- def get_coordinates(self, non_linear=True, method="all", object_cutoff=0, use_flat=False):
+ def get_coordinates(
+ self, non_linear=True, method="all", object_cutoff=0, use_flat=False
+ ):
"""Extracts pixel coordinates from the segmentation data.
Parameters
@@ -178,7 +178,7 @@ class PyNutil:
method=method,
object_cutoff=object_cutoff,
atlas_volume=self.atlas_volume,
- use_flat=use_flat
+ use_flat=use_flat,
)
self.pixel_points = pixel_points
self.centroids = centroids
@@ -221,7 +221,9 @@ class PyNutil:
current_points = None
print("region areas list")
print(self.region_areas_list)
- for pl,cl,ra in zip(self.points_len, self.centroids_len, self.region_areas_list):
+ 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"):
@@ -231,41 +233,52 @@ class PyNutil:
)
print("current df", current_df)
-# create the df for section report and all report
-# pixel_count_per_region returns a df with idx, pixel count, name and RGB.
-# ra is region area list from
-# merge current_df onto ra (region_areas_list) based on idx column
-#(left means use only keys from left frame, preserve key order)
-
+ # create the df for section report and all report
+ # pixel_count_per_region returns a df with idx, pixel count, name and RGB.
+ # ra is region area list from
+ # merge current_df onto ra (region_areas_list) based on idx column
+ # (left means use only keys from left frame, preserve key order)
+
"""
Merge region areas and object areas onto the atlas label file.
Remove duplicate columns
Calculate and add area_fraction to new column in the df.
"""
-
- all_region_df = self.atlas_labels.merge(ra, on = 'idx', how='left')
- current_df_new = all_region_df.merge(current_df, on= 'idx', how= 'left', suffixes= (None,"_y")).drop(columns=["a","VIS", "MSH", "name_y","r_y","g_y","b_y"])
- current_df_new["area_fraction"] = current_df_new["pixel_count"] / current_df_new["region_area"]
+
+ all_region_df = self.atlas_labels.merge(ra, on="idx", how="left")
+ current_df_new = all_region_df.merge(
+ current_df, on="idx", how="left", suffixes=(None, "_y")
+ ).drop(columns=["a", "VIS", "MSH", "name_y", "r_y", "g_y", "b_y"])
+ current_df_new["area_fraction"] = (
+ current_df_new["pixel_count"] / current_df_new["region_area"]
+ )
current_df_new.fillna(0, inplace=True)
per_section_df.append(current_df_new)
prev_pl += pl
prev_cl += cl
-
-
+
##combine all the slice reports, groupby idx, name, rgb and sum region and object pixels. Remove area_fraction column and recalculate.
- self.label_df = pd.concat(per_section_df).groupby(['idx','name','r','g','b']).sum().reset_index().drop(columns=['area_fraction'])
- self.label_df["area_fraction"] = self.label_df["pixel_count"] / self.label_df["region_area"]
+ self.label_df = (
+ pd.concat(per_section_df)
+ .groupby(["idx", "name", "r", "g", "b"])
+ .sum()
+ .reset_index()
+ .drop(columns=["area_fraction"])
+ )
+ self.label_df["area_fraction"] = (
+ self.label_df["pixel_count"] / self.label_df["region_area"]
+ )
self.label_df.fillna(0, inplace=True)
"""
Potential source of error:
If there are duplicates in the label file, regional results will be duplicated and summed leading to incorrect results
"""
- #reorder the df to match the order of idx column in self.atlas_labels
- self.label_df = self.label_df.set_index('idx')
- self.label_df = self.label_df.reindex(index=self.atlas_labels['idx'])
+ # reorder the df to match the order of idx column in self.atlas_labels
+ self.label_df = self.label_df.set_index("idx")
+ self.label_df = self.label_df.reindex(index=self.atlas_labels["idx"])
self.label_df = self.label_df.reset_index()
-
+
self.labeled_points = labeled_points
self.labeled_points_centroids = labeled_points_centroids
self.per_section_df = per_section_df
@@ -289,7 +302,7 @@ class PyNutil:
"""
if not os.path.exists(output_folder):
os.makedirs(output_folder)
-
+
if not os.path.exists(f"{output_folder}/whole_series_report"):
os.makedirs(f"{output_folder}/whole_series_report")
@@ -299,9 +312,11 @@ class PyNutil:
"if you want to save the quantification please run quantify_coordinates"
)
else:
-
self.label_df.to_csv(
- f"{output_folder}/whole_series_report/counts.csv", sep=";", na_rep="", index=False
+ f"{output_folder}/whole_series_report/counts.csv",
+ sep=";",
+ na_rep="",
+ index=False,
)
if not os.path.exists(f"{output_folder}/per_section_meshview"):
os.makedirs(f"{output_folder}/per_section_meshview")
diff --git a/PyNutil/propagation.py b/PyNutil/propagation.py
index d4af1e8..1f3b643 100644
--- a/PyNutil/propagation.py
+++ b/PyNutil/propagation.py
@@ -1,121 +1,132 @@
-import math,re
+import math, re
+
def propagate(arr):
- for slice in arr:
- if "nr" not in slice:
- slice["nr"]=int(re.search(r"_s(\d+)",slice["filename"]).group(1))
-
- arr.sort(key=lambda slice:slice["nr"])
-
- linregs=[LinReg() for i in range(11)]
- count=0
- for slice in arr:
- if "anchoring" in slice:
- a=slice["anchoring"]
- for i in range(3):
- a[i]+=(a[i+3]+a[i+6])/2
- a.extend([normalize(a,3)/slice["width"],normalize(a,6)/slice["height"]])
- for i in range(len(linregs)):
- linregs[i].add(slice["nr"],a[i])
- count+=1
-
- if count>=2:
- l=len(arr)
- if not "anchoring" in arr[0]:
- nr=arr[0]["nr"]
- a=[linreg.get(nr) for linreg in linregs]
- orthonormalize(a)
- arr[0]["anchoring"]=a
- count+=1
- if not "anchoring" in arr[l-1]:
- nr=arr[l-1]["nr"]
- a=[linreg.get(nr) for linreg in linregs]
- orthonormalize(a)
- arr[l-1]["anchoring"]=a
- count+=1
-
- start=1
- while count<l:
- while "anchoring" in arr[start]:
- start+=1
- next=start+1
- while not "anchoring" in arr[next]:
- next+=1
- pnr=arr[start-1]["nr"]
- nnr=arr[next]["nr"]
- panch=arr[start-1]["anchoring"]
- nanch=arr[next]["anchoring"]
- linints=[LinInt(pnr,panch[i],nnr,nanch[i]) for i in range(len(panch))]
- for i in range(start,next):
- nr=arr[i]["nr"]
- arr[i]["anchoring"]=[linint.get(nr) for linint in linints]
- count+=1
- start=next+1
+ for slice in arr:
+ if "nr" not in slice:
+ slice["nr"] = int(re.search(r"_s(\d+)", slice["filename"]).group(1))
+ arr.sort(key=lambda slice: slice["nr"])
+
+ linregs = [LinReg() for i in range(11)]
+ count = 0
for slice in arr:
- a=slice["anchoring"]
- orthonormalize(a)
- v=a.pop()
- u=a.pop()
- for i in range(3):
- a[i+3]*=u*slice["width"]
- a[i+6]*=v*slice["height"]
- a[i]-=(a[i+3]+a[i+6])/2
- return arr
+ if "anchoring" in slice:
+ a = slice["anchoring"].copy()
+ for i in range(3):
+ a[i] += (a[i + 3] + a[i + 6]) / 2
+ a.extend(
+ [normalize(a, 3) / slice["width"], normalize(a, 6) / slice["height"]]
+ )
+ for i in range(len(linregs)):
+ linregs[i].add(slice["nr"], a[i])
+ count += 1
+
+ if count >= 2:
+ l = len(arr)
+ if not "anchoring" in arr[0]:
+ nr = arr[0]["nr"]
+ a = [linreg.get(nr) for linreg in linregs]
+ orthonormalize(a)
+ arr[0]["anchoring"] = a
+ count += 1
+ if not "anchoring" in arr[l - 1]:
+ nr = arr[l - 1]["nr"]
+ a = [linreg.get(nr) for linreg in linregs]
+ orthonormalize(a)
+ arr[l - 1]["anchoring"] = a
+ count += 1
+
+ start = 1
+ while count < l:
+ while "anchoring" in arr[start]:
+ start += 1
+ next = start + 1
+ while not "anchoring" in arr[next]:
+ next += 1
+ pnr = arr[start - 1]["nr"]
+ nnr = arr[next]["nr"]
+ panch = arr[start - 1]["anchoring"]
+ nanch = arr[next]["anchoring"]
+ linints = [LinInt(pnr, panch[i], nnr, nanch[i]) for i in range(len(panch))]
+ for i in range(start, next):
+ nr = arr[i]["nr"]
+ arr[i]["anchoring"] = [linint.get(nr) for linint in linints]
+ count += 1
+ start = next + 1
+
+ for slice in arr:
+ a = slice["anchoring"]
+ orthonormalize(a)
+ v = a.pop()
+ u = a.pop()
+ for i in range(3):
+ a[i + 3] *= u * slice["width"]
+ a[i + 6] *= v * slice["height"]
+ a[i] -= (a[i + 3] + a[i + 6]) / 2
+ return arr
+
class LinInt:
- def __init__(self,x1,y1,x2,y2):
- self.x1=x1
- self.x2=x2
- self.y1=y1
- self.y2=y2
+ def __init__(self, x1, y1, x2, y2):
+ self.x1 = x1
+ self.x2 = x2
+ self.y1 = y1
+ self.y2 = y2
+
+ def get(self, x):
+ return self.y1 + (self.y2 - self.y1) * (x - self.x1) / (self.x2 - self.x1)
- def get(self,x):
- return self.y1+(self.y2-self.y1)*(x-self.x1)/(self.x2-self.x1)
class LinReg:
- def __init__(self):
- self.n=0
- self.Sx=0
- self.Sy=0
- self.Sxx=0
- self.Sxy=0
-
- def add(self,x,y):
- self.n+=1
- self.Sx+=x
- self.Sy+=y
- self.Sxx+=x*x
- self.Sxy+=x*y
- if self.n>=2:
- self.b=(self.n*self.Sxy-self.Sx*self.Sy)/(self.n*self.Sxx-self.Sx*self.Sx)
- self.a=self.Sy/self.n-self.b*self.Sx/self.n
-
- def get(self,x):
- return self.a+self.b*x
-
-def normalize(arr,idx):
- l=0
- for i in range(3):
- l+=arr[idx+i]*arr[idx+i]
- l=math.sqrt(l)
- for i in range(3):
- arr[idx+i]/=l
- return l
+ def __init__(self):
+ self.n = 0
+ self.Sx = 0
+ self.Sy = 0
+ self.Sxx = 0
+ self.Sxy = 0
+
+ def add(self, x, y):
+ self.n += 1
+ self.Sx += x
+ self.Sy += y
+ self.Sxx += x * x
+ self.Sxy += x * y
+ if self.n >= 2:
+ self.b = (self.n * self.Sxy - self.Sx * self.Sy) / (
+ self.n * self.Sxx - self.Sx * self.Sx
+ )
+ self.a = self.Sy / self.n - self.b * self.Sx / self.n
+
+ def get(self, x):
+ return self.a + self.b * x
+
+
+def normalize(arr, idx):
+ l = 0
+ for i in range(3):
+ l += arr[idx + i] * arr[idx + i]
+ l = math.sqrt(l)
+ for i in range(3):
+ arr[idx + i] /= l
+ return l
+
def orthonormalize(arr):
- normalize(arr,3)
- dot=0
- for i in range(3):
- dot+=arr[i+3]*arr[i+6]
- for i in range(3):
- arr[i+6]-=arr[i+3]*dot
- normalize(arr,6)
-
-if __name__=="__main__":
- import json,sys
- with open(sys.argv[1]) as f:
- series=json.load(f)
- propagate(series["slices"])
- with open(sys.argv[2],"w") as f:
- json.dump(series,f)
+ normalize(arr, 3)
+ dot = 0
+ for i in range(3):
+ dot += arr[i + 3] * arr[i + 6]
+ for i in range(3):
+ arr[i + 6] -= arr[i + 3] * dot
+ normalize(arr, 6)
+
+
+if __name__ == "__main__":
+ import json, sys
+
+ with open(sys.argv[1]) as f:
+ series = json.load(f)
+ propagate(series["slices"])
+ with open(sys.argv[2], "w") as f:
+ json.dump(series, f)
diff --git a/PyNutil/read_and_write.py b/PyNutil/read_and_write.py
index ba5420f..8be7f08 100644
--- a/PyNutil/read_and_write.py
+++ b/PyNutil/read_and_write.py
@@ -8,6 +8,7 @@ import nrrd
import re
from .propagation import propagate
+
# related to read and write
# this function reads a VisuAlign JSON and returns the slices
def load_visualign_json(filename):
@@ -24,19 +25,21 @@ def load_visualign_json(filename):
name = os.path.basename(filename)
lz_compat_file = {
- "name":name,
- "target":vafile["atlas"],
- "target-resolution":[ 456, 528, 320],
- "slices":slices,
+ "name": name,
+ "target": vafile["atlas"],
+ "target-resolution": [456, 528, 320],
+ "slices": slices,
}
- # save with .json extension need to see if i can remove this
- with open(filename.replace(".waln", ".json").replace(".wwrp", ".json"), "w") as f:
+ # save with .json extension need to see if i can remove this
+ with open(
+ filename.replace(".waln", ".json").replace(".wwrp", ".json"), "w"
+ ) as f:
json.dump(lz_compat_file, f, indent=4)
else:
slices = vafile["slices"]
-
- slices = propagate(slices)
+ if len(slices) > 1:
+ slices = propagate(slices)
return slices
@@ -95,24 +98,27 @@ def flat_to_array(file, labelfile):
# It has the shape width times height
data = struct.unpack(">" + ("xBH"[b] * (w * h)), f.read(b * w * h))
elif file.endswith(".seg"):
- with open(file,"rb") as f:
+ with open(file, "rb") as f:
+
def byte():
return f.read(1)[0]
+
def code():
- c=byte()
- if(c<0):
+ c = byte()
+ if c < 0:
raise "!"
- return c if c<128 else (c&127)|(code()<<7)
+ return c if c < 128 else (c & 127) | (code() << 7)
+
if "SegRLEv1" != f.read(8).decode():
raise "Header mismatch"
- atlas=f.read(code()).decode()
+ atlas = f.read(code()).decode()
print(f"Target atlas: {atlas}")
- codes=[code() for x in range(code())]
- w=code()
- h=code()
- data=[]
- while len(data)<w*h:
- data += [codes[byte() if len(codes)<=256 else code()]]*(code()+1)
+ codes = [code() for x in range(code())]
+ w = code()
+ h = code()
+ data = []
+ while len(data) < w * h:
+ data += [codes[byte() if len(codes) <= 256 else code()]] * (code() + 1)
# convert flat file data into an array, previously data was a tuple
imagedata = np.array(data)
@@ -124,13 +130,15 @@ def flat_to_array(file, labelfile):
image[y, x] = imagedata[x + y * w]
image_arr = np.array(image)
- #return image_arr
+ # return image_arr
"""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[
+ coordsx, coordsy
+ ] # assign x,y coords from image_array into values
lbidx = labelfile["idx"].values
allen_id_image = lbidx[values.astype(int)]
return allen_id_image
diff --git a/PyNutil/reconstruct_dzi.py b/PyNutil/reconstruct_dzi.py
index 8a43697..bbdff40 100644
--- a/PyNutil/reconstruct_dzi.py
+++ b/PyNutil/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.
@@ -24,15 +25,21 @@ def reconstruct_dzi(zip_file_path):
# Get the filenames of the highest level tiles
highest_level_files = [
- i for i in zip_file.namelist() if i.endswith(".png") and i.split("/")[-2] == highest_level
+ i
+ for i in zip_file.namelist()
+ if i.endswith(".png") and i.split("/")[-2] == highest_level
]
# Read the DZI file
- dzi_file = zip_file.open([i for i in zip_file.namelist() if i.endswith(".dzi")][0])
+ dzi_file = zip_file.open(
+ [i for i in zip_file.namelist() if i.endswith(".dzi")][0]
+ )
xml = dzi_file.read()
json_data = xmltodict.parse(xml)
tileSize = json_data["Image"]["@TileSize"]
- width, height = int(json_data["Image"]["Size"]["@Width"]), int(json_data["Image"]["Size"]["@Height"])
+ width, height = int(json_data["Image"]["Size"]["@Width"]), int(
+ json_data["Image"]["Size"]["@Height"]
+ )
# Create an empty image
image = np.zeros((height, width, 3), dtype=np.uint8)
@@ -48,4 +55,4 @@ def reconstruct_dzi(zip_file_path):
image[y : y + image_.shape[0], x : x + image_.shape[1], :] = image_
# Save the image
- return image
\ No newline at end of file
+ return image
diff --git a/PyNutil/visualign_deformations.py b/PyNutil/visualign_deformations.py
index 534dfa7..2caef65 100644
--- a/PyNutil/visualign_deformations.py
+++ b/PyNutil/visualign_deformations.py
@@ -2,12 +2,13 @@
import numpy as np
from scipy.spatial import Delaunay
+
def triangulate(w, h, markers):
vertices = [
[-0.1 * w, -0.1 * h, -0.1 * w, -0.1 * h],
- [ 1.1 * w, -0.1 * h, 1.1 * w, -0.1 * h],
- [-0.1 * w, 1.1 * h, -0.1 * w, 1.1 * h],
- [ 1.1 * w, 1.1 * h, 1.1 * w, 1.1 * h],
+ [1.1 * w, -0.1 * h, 1.1 * w, -0.1 * h],
+ [-0.1 * w, 1.1 * h, -0.1 * w, 1.1 * h],
+ [1.1 * w, 1.1 * h, 1.1 * w, 1.1 * h],
]
edges = [0] * ((len(markers) + 4) * (len(markers) + 4 - 1) // 2)
@@ -155,7 +156,7 @@ class Triangle:
[[bx - ax, by - ay, 0], [cx - ax, cy - ay, 0], [ax, ay, 1]]
)
if self.decomp == None:
- print('e')
+ print("e")
a2 = distsquare(bx, by, cx, cy)
b2 = distsquare(ax, ay, cx, cy)
c2 = distsquare(ax, ay, bx, by)
@@ -190,7 +191,6 @@ class Triangle:
# xy: 2-dimensional array with one xy-pair per row
def inforward_vec(self, x, y, xPrime, yPrime):
-
uv1 = rowmul3_vec(x, y, self.forwarddecomp)
# also compute the next step, since it uses the parameters of this triangle
ok = (
@@ -213,7 +213,7 @@ class Triangle:
def intriangle_vec(self, x, y, xPrime, yPrime):
if self.decomp is None:
- print('e')
+ print("e")
uv1 = rowmul3_vec(x, y, self.decomp)
# also compute the next step, since it uses the parameters of this triangle
ok = (
@@ -233,5 +233,3 @@ class Triangle:
+ (self.B[1] - self.A[1]) * uv1[ok, 0]
+ (self.C[1] - self.A[1]) * uv1[ok, 1]
)
-
-
diff --git a/messing_around_files/pyflat.py b/messing_around_files/pyflat.py
index cbe6605..f4fbfb5 100644
--- a/messing_around_files/pyflat.py
+++ b/messing_around_files/pyflat.py
@@ -1,17 +1,17 @@
-
"""This is Gergely Csusc's code from VisuAlign MediaWiki NITRC page,
-this was a test, not used elsewhere in PyNutil"""
+this was a test, not used elsewhere in PyNutil"""
import sys
+
for arg in sys.argv:
if arg.startswith("flat="):
- flatfile=arg[len("flat="):]
+ flatfile = arg[len("flat=") :]
if arg.startswith("label="):
- labelfile=arg[len("label="):]
+ labelfile = arg[len("label=") :]
if arg.startswith("json="):
- jsonfile=arg[len("json="):]
+ jsonfile = arg[len("json=") :]
if arg.startswith("output="):
- outfile=arg[len("output="):]
+ outfile = arg[len("output=") :]
if "flatfile" not in vars():
print("flat=<some .flat file> parameter missing")
sys.exit()
@@ -19,31 +19,37 @@ if "outfile" not in vars():
print("output=<new .png file> parameter missing")
sys.exit()
-palette=False
+palette = False
if "labelfile" in vars():
import re
- palette=[]
+
+ palette = []
with open(labelfile) as f:
for line in f:
- lbl=re.match(r'\s*\d+\s+(\d+)\s+(\d+)\s+(\d+)',line)
+ lbl = re.match(r"\s*\d+\s+(\d+)\s+(\d+)\s+(\d+)", line)
if lbl:
- palette.append((int(lbl[1]),int(lbl[2]),int(lbl[3])))
+ palette.append((int(lbl[1]), int(lbl[2]), int(lbl[3])))
print(f"{len(palette)} labels parsed")
elif "jsonfile" in vars():
import json
+
with open(jsonfile) as f:
- palette=[(i["red"],i["green"],i["blue"]) for i in json.load(f)]
+ palette = [(i["red"], i["green"], i["blue"]) for i in json.load(f)]
print(f"{len(palette)} labels loaded")
import struct
-with open(flatfile,"rb") as f:
- b,w,h=struct.unpack(">BII",f.read(9))
- data=struct.unpack(">"+("xBH"[b]*(w*h)),f.read(b*w*h))
+
+with open(flatfile, "rb") as f:
+ b, w, h = struct.unpack(">BII", f.read(9))
+ data = struct.unpack(">" + ("xBH"[b] * (w * h)), f.read(b * w * h))
print(f"{b} bytes per pixel, {w} x {h} resolution")
import PIL.Image
-image=PIL.Image.new("RGB" if palette else "L",(w,h))
+
+image = PIL.Image.new("RGB" if palette else "L", (w, h))
for y in range(h):
for x in range(w):
- image.putpixel((x,y),palette[data[x+y*w]] if palette else data[x+y*w] & 255)
-image.save(outfile,"PNG")
\ No newline at end of file
+ image.putpixel(
+ (x, y), palette[data[x + y * w]] if palette else data[x + y * w] & 255
+ )
+image.save(outfile, "PNG")
diff --git a/messing_around_files/testing_openflatfile.py b/messing_around_files/testing_openflatfile.py
index 1b076ba..a52502b 100644
--- a/messing_around_files/testing_openflatfile.py
+++ b/messing_around_files/testing_openflatfile.py
@@ -39,15 +39,15 @@ with open(base, "rb") as f:
"""assign label file values into image array"""
labelfile = pd.read_csv(r"annotation_volumes\allen2017_colours.csv")
-#labelfile[]
-print(list(zip(val,count)))
+# labelfile[]
+print(list(zip(val, count)))
-temp_copy = labelfile.loc[val.astype(int),:].copy()
+temp_copy = labelfile.loc[val.astype(int), :].copy()
-temp_copy['region areas'] = count
+temp_copy["region areas"] = count
-temp_copy.to_csv('../outputs/region_areas_quick_test.csv',sep=';')
+temp_copy.to_csv("../outputs/region_areas_quick_test.csv", sep=";")
allen_id_image = np.zeros((h, w)) # create an empty image array
diff --git a/scripts/create_test_data.py b/scripts/create_test_data.py
index 1e196fe..4746729 100644
--- a/scripts/create_test_data.py
+++ b/scripts/create_test_data.py
@@ -3,15 +3,16 @@ from PIL import Image, ImageDraw
def generate_test_image(width, height, file_path):
# Create an empty white image
- image = Image.new(mode = "RGB", size = (width, height), color=(255,255,255))
+ image = Image.new(mode="RGB", size=(width, height), color=(255, 255, 255))
# Save the image as a PNG file
image.save(file_path)
+
# Example usage
width = 1000 # Specify the width of the image in pixels
height = 1000 # Specify the height of the image in pixels
-file_path = 'testimage.png' # Specify the file path for saving the image
+file_path = "testimage.png" # Specify the file path for saving the image
# Generate and save the black and white image
generate_test_image(width, height, file_path)
@@ -19,10 +20,13 @@ generate_test_image(width, height, file_path)
"""This is used to generate the test data"""
-def generate_image_with_squares(width, height, square_diameter, square_locations, num_images):
+
+def generate_image_with_squares(
+ width, height, square_diameter, square_locations, num_images
+):
# Create a white image with objects of specified size at specified x, y locations
for i in range(1, num_images + 1):
- image = Image.new('RGB', (width, height), color=(255, 255, 255))
+ image = Image.new("RGB", (width, height), color=(255, 255, 255))
# Create a draw object
draw = ImageDraw.Draw(image)
@@ -37,13 +41,22 @@ def generate_image_with_squares(width, height, square_diameter, square_locations
file_name = f"../test_data/PyTest/test_s00{i}.png"
image.save(file_name, "PNG")
+
# Example usage
width = 1500 # Specify the width of the image in pixels
-height = 1000 # Specify the height of the image in pixels
+height = 1000 # Specify the height of the image in pixels
square_diameter = 10 # Specify the size of the black squares in pixels
-square_locations = [(500, 500), (500, 600), (500, 700), (1000, 500), (1000,600),(1000,700)] # Specify the x, y locations of the black squares
+square_locations = [
+ (500, 500),
+ (500, 600),
+ (500, 700),
+ (1000, 500),
+ (1000, 600),
+ (1000, 700),
+] # Specify the x, y locations of the black squares
num_images = 5
# Generate the white image with black squares
-image = generate_image_with_squares(width, height, square_diameter, square_locations,num_images)
-
+image = generate_image_with_squares(
+ width, height, square_diameter, square_locations, num_images
+)
diff --git a/testOOP.py b/testOOP.py
index 8075fb8..368f2c4 100644
--- a/testOOP.py
+++ b/testOOP.py
@@ -1,7 +1,7 @@
from PyNutil import PyNutil
pnt = PyNutil(settings_file=r"PyNutil/test/test8_PyNutil_fixed.json")
-##Use flat can be set to True if you want to use the flat file
+##Use flat can be set to True if you want to use the flat file
# instead of the visualign json (this is only useful for testing and will be removed)
pnt.get_coordinates(object_cutoff=0, use_flat=False)
@@ -10,4 +10,4 @@ pnt.quantify_coordinates()
pnt.save_analysis("PyNutil/outputs/test8_PyNutil")
# remove name, r, g, b, from pixel_
-# add to region_areas df
\ No newline at end of file
+# add to region_areas df
--
GitLab