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/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 ab1365d4d7ce2d3a5cf274d8be9b2ffd120833a4..c6369444ea40279ad1a26b71e25f5ccfc639449e 100644
--- a/testOOP.py
+++ b/testOOP.py
@@ -1,9 +1,11 @@
from PyNutil import PyNutil
+import os
+os.chdir("..")
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
# instead of the visualign json (this is only useful for testing and will be removed)
-pnt.get_coordinates(object_cutoff=0, use_flat=True)
+pnt.get_coordinates(object_cutoff=0, use_flat=False)
pnt.quantify_coordinates()