diff --git a/.gitignore b/.gitignore
index 67f7b335606483f98974953b31db180a357277d1..7ad8120dce388238b14d68eb520b94aa1977e489 100644
--- a/.gitignore
+++ b/.gitignore
@@ -164,5 +164,6 @@ cython_debug/
test_data/tTA_2877_NOP
test_data/tTA_2877_NOP_horizontal_final_2017.json
+fileForMessingAround.py
-fileForMessingAround.py
\ No newline at end of file
+outputs/
diff --git a/PyNutil/coordinate_extraction.py b/PyNutil/coordinate_extraction.py
index 73a326a0441e07dff78c9da4cf5313b9e5cea69b..a58ae144afc2d95f0050d21a168a0b2c7af237f2 100644
--- a/PyNutil/coordinate_extraction.py
+++ b/PyNutil/coordinate_extraction.py
@@ -234,9 +234,55 @@ def folder_to_atlas_space(
segmentations,
)
+def load_segmentation(segmentation_path: str):
+ """Load a segmentation from a file."""
+ print(f"working on {segmentation_path}")
+ if segmentation_path.endswith(".dzip"):
+ print("Reconstructing dzi")
+ return reconstruct_dzi(segmentation_path)
+ else:
+ return cv2.imread(segmentation_path)
+
+def remove_background(segmentation: np.array):
+ """Remove the background from the segmentation and return the pixel id."""
+ segmentation_no_background = segmentation[~np.all(segmentation == 0, axis=2)]
+ print("length of non background pixels: ", len(segmentation_no_background))
+ pixel_id = segmentation_no_background[0]
+ print("detected pixel_id: ", pixel_id)
+ return segmentation_no_background, pixel_id
+
+def get_region_areas(use_flat, atlas_labels, flat_file_atlas, seg_width, seg_height, slice_dict, atlas_volume, triangulation):
+ if use_flat:
+ 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_dict["anchoring"],
+ atlas_volume,
+ triangulation
+ )
+ return region_areas
+
+def get_transformed_coordinates(non_linear, slice_dict, method, scaled_x, scaled_y, centroids, scaled_centroidsX, scaled_centroidsY, triangulation):
+ new_x, new_y, centroids_new_x, centroids_new_y = None, None, None, None
+ if non_linear and "markers" in slice_dict:
+ if method in ["per_pixel", "all"] and scaled_x is not None:
+ new_x, new_y = transform_vec(triangulation, scaled_x, scaled_y)
+ if method in ["per_object", "all"] and centroids is not None:
+ centroids_new_x, centroids_new_y = transform_vec(triangulation, scaled_centroidsX, scaled_centroidsY)
+ else:
+ if method in ["per_pixel", "all"]:
+ new_x, new_y = scaled_x, scaled_y
+ if method in ["per_object", "all"]:
+ centroids_new_x, centroids_new_y = scaled_centroidsX, scaled_centroidsY
+ return new_x, new_y, centroids_new_x, centroids_new_y
def segmentation_to_atlas_space(
- slice,
+ slice_dict,
segmentation_path,
atlas_labels,
flat_file_atlas=None,
@@ -251,113 +297,30 @@ def segmentation_to_atlas_space(
atlas_volume=None,
use_flat=False,
):
- """Combines many functions to convert a segmentation to atlas space. It takes care
- of deformations."""
- print(f"working on {segmentation_path}")
- if segmentation_path.endswith(".dzip"):
- print("Reconstructing dzi")
- segmentation = reconstruct_dzi(segmentation_path)
-
- else:
- segmentation = cv2.imread(segmentation_path)
+ segmentation = load_segmentation(segmentation_path)
if pixel_id == "auto":
-
- # Remove the background from the segmentation
- segmentation_no_background = segmentation[~np.all(segmentation == 0, axis=2)]
- # pixel_id = np.vstack(
- # {tuple(r) for r in segmentation_no_background.reshape(-1, 3)}
- # ) # Remove background
- # 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)
-
- # Transform pixels to registration space (the registered image and segmentation have different dimensions)
- seg_height = segmentation.shape[0]
- seg_width = segmentation.shape[1]
- 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)
- )
+ segmentation, pixel_id = remove_background(segmentation)
+ seg_height, seg_width = segmentation.shape[:2]
+ reg_height, reg_width = slice_dict["height"], slice_dict["width"]
+ if non_linear and "markers" in slice_dict:
+ triangulation = triangulate(reg_width, reg_height, slice_dict["markers"])
else:
- 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
- )
-
+ triangulation = None
+ region_areas = get_region_areas(use_flat, atlas_labels, flat_file_atlas, seg_width, seg_height, slice_dict, atlas_volume, triangulation)
+ y_scale, x_scale = transform_to_registration(seg_height, seg_width, reg_height, reg_width)
centroids, points = None, None
-
if method in ["per_object", "all"]:
- centroids, scaled_centroidsX, scaled_centroidsY = get_centroids(
- segmentation, pixel_id, y_scale, x_scale, object_cutoff
- )
+ centroids, scaled_centroidsX, scaled_centroidsY = get_centroids(segmentation, pixel_id, y_scale, x_scale, object_cutoff)
if method in ["per_pixel", "all"]:
scaled_y, scaled_x = get_scaled_pixels(segmentation, pixel_id, y_scale, x_scale)
- 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"]:
- if scaled_x is not None:
- new_x, new_y = transform_vec(triangulation, scaled_x, scaled_y)
- else:
- new_x, new_y = scaled_x, scaled_y
- if method in ["per_object", "all"]:
- if centroids is not None:
- centroids_new_x, centroids_new_y = transform_vec(
- triangulation, scaled_centroidsX, scaled_centroidsY
- )
- else:
- 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."
- )
- if method in ["per_pixel", "all"]:
- new_x, new_y = scaled_x, scaled_y
- if method in ["per_object", "all"]:
- centroids_new_x, centroids_new_y = scaled_centroidsX, scaled_centroidsY
- else:
- if method in ["per_pixel", "all"]:
- new_x, new_y = scaled_x, scaled_y
- if method in ["per_object", "all"]:
- centroids_new_x, centroids_new_y = scaled_centroidsX, scaled_centroidsY
- # Scale U by Uxyz/RegWidth and V by Vxyz/RegHeight
- if method in ["per_pixel", "all"]:
- if new_x is not None:
- points = transform_to_atlas_space(
- slice["anchoring"], new_y, new_x, reg_height, reg_width
- )
- else:
- points = np.array([])
- 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,
- )
- else:
- centroids = np.array([])
-
-
- points_list[index] = np.array(points)
- centroids_list[index] = np.array(centroids)
+ new_x, new_y, centroids_new_x, centroids_new_y = get_transformed_coordinates(non_linear, slice_dict, method, scaled_x, scaled_y, centroids, scaled_centroidsX, scaled_centroidsY, triangulation)
+ if method in ["per_pixel", "all"] and new_x is not None:
+ points = transform_to_atlas_space(slice_dict["anchoring"], new_y, new_x, reg_height, reg_width)
+ if method in ["per_object", "all"] and centroids_new_x is not None:
+ centroids = transform_to_atlas_space(slice_dict["anchoring"], centroids_new_y, centroids_new_x, reg_height, reg_width)
+ points_list[index] = np.array(points if points is not None else [])
+ centroids_list[index] = np.array(centroids if centroids is not None else [])
region_areas_list[index] = region_areas
diff --git a/PyNutil/counting_and_load.py b/PyNutil/counting_and_load.py
index 94c151a70edfb5dc40d62dc87b1f3284f55f0aa3..50baba96ea65d8d0be849b29a3468c79746d2cbb 100644
--- a/PyNutil/counting_and_load.py
+++ b/PyNutil/counting_and_load.py
@@ -3,7 +3,7 @@ import pandas as pd
import struct
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):
@@ -117,10 +117,7 @@ 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"""
-import struct
-import cv2
-import numpy as np
-import pandas as pd
+
def read_flat_file(file):
@@ -187,12 +184,41 @@ def count_pixels_per_label(image, scale_factor=False):
return df_area_per_label
+def warp_image(image, triangulation, rescaleXY):
+ if rescaleXY is not None:
+ w,h = rescaleXY
+ else:
+ h, w = image.shape
+ reg_h, reg_w = image.shape
+ oldX, oldY = np.meshgrid(np.arange(reg_w), np.arange(reg_h))
+ oldX = oldX.flatten()
+ oldY = oldY.flatten()
+ h_scale = h / reg_h
+ w_scale = w / reg_w
+ oldX = oldX * w_scale
+ oldY = oldY * h_scale
+ newX, newY = transform_vec(triangulation, oldX, oldY)
+ newX = newX / w_scale
+ newY = newY / h_scale
+ newX = newX.reshape(reg_h, reg_w)
+ newY = newY.reshape(reg_h, reg_w)
+ newX = newX.astype(int)
+ newY = newY.astype(int)
+ newX[newX >= reg_w] = reg_w - 1
+ newY[newY >= reg_h] = reg_h - 1
+ newX[newX < 0] = 0
+ newY[newY < 0] = 0
+ new_image = image[newY, newX]
+ return new_image
+
def flat_to_dataframe(
- labelfile, file=None, rescaleXY=None, image_vector=None, volume=None
+ labelfile, file=None, rescaleXY=None, image_vector=None, volume=None, triangulation=None
):
if (image_vector is not None) and (volume is not None):
image = generate_target_slice(image_vector, volume)
image = np.float64(image)
+ if triangulation is not None:
+ image = warp_image(image, triangulation, rescaleXY)
elif file.endswith(".flat"):
image = read_flat_file(file)
elif file.endswith(".seg"):
diff --git a/PyNutil/generate_target_slice.py b/PyNutil/generate_target_slice.py
index 31681b869aee92ba1943a1890cb7bb74556915e8..f042f707bb39e8bca8bbdef845af0af3c8ffb6c5 100644
--- a/PyNutil/generate_target_slice.py
+++ b/PyNutil/generate_target_slice.py
@@ -1,83 +1,35 @@
import numpy as np
+import math
+def generate_target_slice(ouv, atlas):
+ width = None
+ height = None
+ ox, oy, oz, ux, uy, uz, vx, vy, vz = ouv
+ width = np.floor(math.hypot(ux,uy,uz)).astype(int) + 1
+ height = np.floor(math.hypot(vx,vy,vz)).astype(int) + 1
+ data = np.zeros((width, height), dtype=np.uint32).flatten()
+ xdim, ydim, zdim = atlas.shape
+ y_values = np.arange(height)
+ x_values = np.arange(width)
+ hx = ox + vx * (y_values / height)
+ hy = oy + vy * (y_values / height)
+ hz = oz + vz * (y_values / height)
+ wx = ux * (x_values / width)
+ wy = uy * (x_values / width)
+ wz = uz * (x_values / width)
+ lx = np.floor(hx[:, None] + wx).astype(int)
+ ly = np.floor(hy[:, None] + wy).astype(int)
+ lz = np.floor(hz[:, None] + wz).astype(int)
+ valid_indices = (0 <= lx) & (lx < xdim) & (0 <= ly) & (ly < ydim) & (0 <= lz) & (lz < zdim)
+ valid_indices = valid_indices.flatten()
+ lxf = lx.flatten()
+ lyf = ly.flatten()
+ lzf = lz.flatten()
+ valid_lx = lxf[valid_indices]
+ valid_ly = lyf[valid_indices]
+ valid_lz = lzf[valid_indices]
+ atlas_slice = atlas[valid_lx,valid_ly,valid_lz]
+ data[valid_indices] = atlas_slice
+ data_im = data.reshape((height, width))
+ return data_im
-def generate_target_slice(alignment, volume):
- Ox, Oy, Oz, Ux, Uy, Uz, Vx, Vy, Vz = alignment
- ##just for mouse for now
- bounds = [455, 527, 319]
- X_size = np.sqrt(np.sum(np.square((Ux, Uy, Uz))))
- Z_size = np.sqrt(np.sum(np.square((Vx, Vy, Vz))))
- X_size = np.round(X_size).astype(int)
- Z_size = np.round(Z_size).astype(int)
- # make this into a grid (0,0) to (320,456)
- Uarange = np.arange(0, 1, 1 / X_size)
- Varange = np.arange(0, 1, 1 / Z_size)
- Ugrid, Vgrid = np.meshgrid(Uarange, Varange)
- Ugrid_x = Ugrid * Ux
- Ugrid_y = Ugrid * Uy
- Ugrid_z = Ugrid * Uz
- Vgrid_x = Vgrid * Vx
- Vgrid_y = Vgrid * Vy
- Vgrid_z = Vgrid * Vz
-
- X_Coords = (Ugrid_x + Vgrid_x).flatten() + Ox
- Y_Coords = (Ugrid_y + Vgrid_y).flatten() + Oy
- Z_Coords = (Ugrid_z + Vgrid_z).flatten() + Oz
-
- X_Coords = np.round(X_Coords).astype(int)
- Y_Coords = np.round(Y_Coords).astype(int)
- 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_pad = X_Coords.copy()
- Y_pad = Y_Coords.copy()
- Z_pad = Z_Coords.copy()
-
- X_pad[out_bounds_Coords] = 0
- Y_pad[out_bounds_Coords] = 0
- Z_pad[out_bounds_Coords] = 0
-
- regions = volume[X_pad, Y_pad, Z_pad]
- ##this is a quick hack to solve rounding errors
- C = len(regions)
- compare = C - X_size * Z_size
- if abs(compare) == X_size:
- if compare > 0:
- Z_size += 1
- if compare < 0:
- Z_size -= 1
- elif abs(C - X_size * Z_size) == Z_size:
- if compare > 0:
- X_size += 1
- if compare < 0:
- X_size -= 1
- elif abs(C - X_size * Z_size) == Z_size + X_size:
- if compare > 0:
- X_size += 1
- Z_size += 1
- if compare < 0:
- X_size -= 1
- Z_size -= 1
- elif abs(C - X_size * Z_size) == Z_size - X_size:
- if compare > 0:
- X_size += 1
- Z_size -= 1
- if compare < 0:
- X_size -= 1
- Z_size += 1
- elif abs(C - X_size * Z_size) == X_size - Z_size:
- if compare > 0:
- X_size -= 1
- Z_size += 1
- if compare < 0:
- X_size += 1
- Z_size -= 1
- regions = regions.reshape((abs(Z_size), abs(X_size)))
- return regions
diff --git a/README.md b/README.md
index e46eb9168fef8f01622def51d3a0680b98d92899..14db1651be2571ec0e8a0425b70485f086edccb2 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,30 @@
# PyNutil
-PyNutil is a Python library for brain-wide quantification and spatial analysis of features in histological images from mouse and rat brain. It aims to replicate the Quantifier feature of the Nutil software (RRID: SCR_017183), to be intergrated with the other QUINT workflow web-tools. It builds on registration to a standardised reference atlas with the QuickNII (RRID:SCR_016854) and VisuAlign software (RRID:SCR_017978) and feature extraction by segmentation with an image analysis software such as ilastik (RRID:SCR_015246).
+PyNutil is a Python library for brain-wide quantification and spatial analysis of features in serial section images from mouse and rat brain. It aims to replicate the Quantifier feature of the Nutil software (RRID: SCR_017183) to be run locally or to be intergrated with the QUINT workflow web-tools. It builds on registration to a standardised reference atlas with the QuickNII (RRID:SCR_016854) and VisuAlign software (RRID:SCR_017978) and feature extraction by segmentation with an image analysis software such as ilastik (RRID:SCR_015246).
For more information about the QUINT workflow:
https://quint-workflow.readthedocs.io/en/latest/
+
+# Usage
+As input PyNutil requires:
+1. An alignment JSON generated with the QuickNII or VisuAlign software
+2. A segmentation file for each brain section with the feature-of-interests displayed in a unique RGB colour (it currently accepts many image formats, png, jpg, jpeg, etc).
+
+To run PyNutil, first fill in a test.json with the correct reference atlas and paths to the required input files. e.g. reference atlas volume, atlas label file, segmentation directory, and path to the alignment json (from QuickNII or VisuAlign).
+
+Then Run testOOP.py outside of the PyNutil directory (cd ..) to inititate the job.
+
+```
+from PyNutil import PyNutil
+
+pnt = PyNutil(settings_file=r"PyNutil/test/test.json")
+
+pnt.get_coordinates(object_cutoff=0, use_flat=False)
+
+pnt.quantify_coordinates()
+
+pnt.save_analysis("PyNutil/outputs/myResults")
+```
+PyNutil generates a series of reports which are saved to the "outputs" directory.
# Acknowledgements
PyNutil is developed at the Neural Systems Laboratory at the Institute of Basic Medical Sciences, University of Oslo, Norway with support from the EBRAINS infrastructure, and funding support from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Framework Partnership Agreement No. 650003 (HBP FPA).
diff --git a/messing_around_files/demo_image_warp.py b/messing_around_files/demo_image_warp.py
new file mode 100644
index 0000000000000000000000000000000000000000..17d55c51fde4aa296a6e84e35c179f3ac3cbdf23
--- /dev/null
+++ b/messing_around_files/demo_image_warp.py
@@ -0,0 +1,67 @@
+import json
+import cv2
+import matplotlib.pyplot as plt
+import os
+import sys
+import numpy as np
+import nrrd
+atlas_path = r"/home/harryc/Github/PyNutilWeb/server/PyNutil/PyNutil/metadata/annotation_volumes/annotation_25_reoriented_2017.nrrd"
+
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+from PyNutil.generate_target_slice import generate_target_slice
+from PyNutil.visualign_deformations import triangulate, transform_vec
+
+def make_slice_ordinal(data):
+ for i, slice in enumerate(np.unique(data)):
+ data[data==slice] = i
+ return data
+data_path = r"/home/harryc/Github/PyNutilWeb/server/PyNutil/test_data/PyNutil_testdataset_Nonlin_SY_fixed_bigcaudoputamen.json"
+with open(data_path, "r") as f:
+ data = json.load(f)
+
+volume, _ = nrrd.read(atlas_path)
+demo = data["slices"][0]
+demo_alignment = demo["anchoring"]
+demo_markers = demo["markers"]
+h = demo["height"]
+w = demo["width"]
+
+
+image = generate_target_slice(demo_alignment, volume)
+image = make_slice_ordinal(image)
+plt.imshow(image)
+plt.show()
+
+triangulation = triangulate(w, h, demo_markers)
+
+
+
+
+
+def warp_image(image, triangulation, h,w):
+ reg_h, reg_w = image.shape
+
+ oldX, oldY = np.meshgrid(np.arange(reg_w), np.arange(reg_h))
+ oldX = oldX.flatten()
+ oldY = oldY.flatten()
+ h_scale = h / reg_h
+ w_scale = w / reg_w
+ oldX = oldX * w_scale
+ oldY = oldY * h_scale
+ newX, newY = transform_vec(triangulation, oldX, oldY)
+ newX = newX / w_scale
+ newY = newY / h_scale
+ newX = newX.reshape(reg_h, reg_w)
+ newY = newY.reshape(reg_h, reg_w)
+ newX = newX.astype(int)
+ newY = newY.astype(int)
+ newX[newX >= reg_w] = reg_w - 1
+ newY[newY >= reg_h] = reg_h - 1
+ newX[newX < 0] = 0
+ newY[newY < 0] = 0
+ new_image = image[newY, newX]
+ return new_image
+
+new_image = warp_image(image, triangulation, h, w)
+plt.imshow(new_image)
+plt.show()
\ No newline at end of file
diff --git a/testOOP.py b/testOOP.py
index 236c95192ee0fe24b243aed406b10e91062c2523..ecf26b3ee19822fc651d80e369b3a56eefa409a3 100644
--- a/testOOP.py
+++ b/testOOP.py
@@ -1,6 +1,9 @@
from PyNutil import PyNutil
+import os
+os.chdir("..")
pnt = PyNutil(settings_file=r"PyNutil/test/test9_PyNutil_linear_only.json")
+
##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)