completed the workflow and added a label for outside the atlas

PyNutil.py

@@ -95,7 +95,7 @@ def SegmentationToAtlasSpace(slice, SegmentationPath, pixelID='auto', nonLinear=
     #scale U by Uxyz/RegWidth and V by Vxyz/RegHeight
     points = transformToAtlasSpace(slice['anchoring'], newY, newX, RegHeight, RegWidth)
     # points = points.reshape(-1)
-    return points
+    return np.array(points)
 
 
 def FolderToAtlasSpace(folder, QUINT_alignment, pixelID=[0, 0, 0], nonLinear=True):
@@ -113,18 +113,22 @@ def FolderToAtlasSpace(folder, QUINT_alignment, pixelID=[0, 0, 0], nonLinear=Tru
         current_slice = slices[current_slice_index[0][0]]
         ##this converts the segmentation to a point cloud
         points.extend(SegmentationToAtlasSpace(current_slice, SegmentationPath, pixelID, nonLinear))
-    return points
+    return np.array(points)
+def createRegionDict(points, regions):
+    """points is a list of points and regions is an id for each point"""
+    regionDict = {region:points[regions==region].flatten().tolist() for region in np.unique(regions)}
+    return regionDict
 
 def WritePoints(pointsDict, filename, infoFile):
     meshview = [
     {
         "idx": idx,
         "count": len(pointsDict[name])//3,
-        "name"  : str(name),
-        "triplets": str(infoFile["name"].values[0]),
-        "r": str(infoFile["r"].values[0]),
-        "g": str(infoFile["g"].values[0]),
-        "b": str(infoFile["b"].values[0])
+        "name"  :str(infoFile["name"].values[infoFile["allenID"]==name][0]),
+        "triplets": pointsDict[name],
+        "r": str(infoFile["r"].values[infoFile["allenID"]==name][0]),
+        "g": str(infoFile["g"].values[infoFile["allenID"]==name][0]),
+        "b": str(infoFile["b"].values[infoFile["allenID"]==name][0])
     }
     for name, idx in zip(pointsDict.keys(), range(len(pointsDict.keys())))
     ]
@@ -132,6 +136,9 @@ def WritePoints(pointsDict, filename, infoFile):
     with open(filename, "w") as f:
         json.dump(meshview, f)
 
+def WritePointsToMeshview(points, pointNames, filename, infoFile):
+    regionDict = createRegionDict(points, pointNames)
+    WritePoints(regionDict, filename, infoFile)
 
 def labelPoints(points, label_volume, scale_factor=1):
     #first convert the points to 3 columns

basic_workflow.py

@@ -3,12 +3,13 @@ import pandas as pd
 #nrrd just lets us open nrrd files
 import nrrd
 #import our function for converting a folder of segmentations to points
-from PyNutil import FolderToAtlasSpace, labelPoints
-
+from PyNutil import FolderToAtlasSpace, labelPoints, WritePointsToMeshview
 
 
+segmentation_folder = "test_data/oneSection15/"
+alignment_json = "test_data/C68_nonlinear.json"
 #now we can use our function to convert the folder of segmentations to points
-points = FolderToAtlasSpace("test_data/oneSection15/", "test_data/C68_nonlinear.json", pixelID=[255, 0, 255], nonLinear=True)
+points = FolderToAtlasSpace(segmentation_folder,alignment_json, pixelID=[255, 0, 255], nonLinear=True)
 #first we need to find the label file for the volume
 label_path = "annotation_volumes//allen2022_colours.csv"
 #then the path to the volume
@@ -18,3 +19,6 @@ label_df = pd.read_csv(label_path)
 #read the annotation volume, it also has a header but we don't need it
 data, header = nrrd.read(volume_path)
 #now we can get the labels for each point
+labels = labelPoints(points, data, scale_factor=2.5)
+#save points to a meshview json
+WritePointsToMeshview(points, labels,"outputs/points.json", label_df)

scripts/reformat_allen_label.py

@@ -7,5 +7,6 @@ path = "itksnap_label_description.txt"
 #set column names
 df = pd.read_csv(path, sep=" ", header=None,  names=["id", "r", "g", "b", "1a", "1b", "1c", "name"])
 df[["name", "allenID"]] = df["name"].str.split(' - ', expand=True)
+df = df.append({"allenID": 0, "name": "background", "r": 255, "g": 255, "b": 255}, ignore_index=True)
 df.to_csv("allen2022_colours.csv", index=False)
 

testMemoVec.py

-import time
-startTime = time.time()
-
-import numpy as np
-from VisuAlignWarpVec import triangulate,  forwardtransform_vec, forwardtransform
-import json
-from tqdm import tqdm
-import matplotlib.pyplot as plt
-import matplotlib
-matplotlib.use('TkAgg')
-
-from ProjectSegmentation3D_vec import FolderToAtlasSpace, getCentroidsAndArea, assignPointsToRegions
-
-import cv2
-
-# Segmentation = cv2.imread("ext-d000033_PVMouseExtraction_pub-Nutil_Quantifier_analysis-81264-Input_dir/ext-d000009_PVMouse_81264_Samp1_resize15__s013_thumbnail_FinalSegmentation.png")
-#count the number of isolated segments
-# Segmentation[~np.all(Segmentation==255, axis=2)] = 0
-
-# centroids, area = getCentroidsAndArea(Segmentation, 4)
-
-    
-
-
-
-#open colour txt file
-path = "itksnap_label_description.txt"
-
-# convert to pandas dataframe
-import pandas as pd
-# use " " as separator
-#set column names
-df = pd.read_csv(path, sep=" ", header=None,  names=["id", "r", "g", "b", "1a", "1b", "1c", "name"])
-df[["name", "allenID"]] = df["name"].str.split(' - ', expand=True)
-df.to_csv("allen2022_colours.csv", index=False)
-
-
-
-
-# `read the annotation volume annotation_10.nrrd`
-import nrrd
-
-data, header = nrrd.read('annotation_10.nrrd')
-
-
-
-
-
-
-points = FolderToAtlasSpace("oneSection15/", "C68_nonlinear.json", pixelID=[255, 0, 255], nonLinear=True)
-
-
-points = np.reshape(points, (-1, 3))
-
-
-Points10um = points*2.5
-
-SwapData = np.transpose(data, (2,0,1))
-SwapData = SwapData[:, ::-1, ::-1]
-Regions = SwapData[Points10um[:,0].astype(int), Points10um[:,1].astype(int), Points10um[:,2].astype(int)]
-#efficiently allocate points to the corresponding region in a dictionary with region as key
-regionDict = {region : [points[Regions==region]] for region in np.unique(Regions)}
-
-#show frequency of regions
-meshview = []
-points_all = []
-idx = 0
-for region in tqdm(regionDict):
-    temp_points = np.array(regionDict[region])
-    temp_points = temp_points.reshape(-1)
-    #write meshview json
-    if region == 0:
-        infoRow = pd.DataFrame([{"allenID": 0, "name": "background", "r": 255, "g": 255, "b": 255}])
-    else:
-        infoRow = df.loc[df['allenID'] == str(region)]
-    points_all.extend(temp_points.tolist())
-    # meshview.append({
-    #     "idx": str(1),
-    #     "count": str(len(temp_points)//3),
-    #     # "name"  : str(infoRow["name"].values[0]),
-    #     "name":"1",
-    #     "triplets": temp_points.tolist(),
-    #     # "r": str(infoRow['r'].values[0]),
-    #     # "g": str(infoRow['g'].values[0]),
-    #     # "b": str(infoRow['b'].values[0])
-    #     "r": 0,
-    #     "g": 0,
-    #     "b": 255
-    # })
-    # idx += 1
-    #write meshview json
-
-
-    # executionTime = (time.time() - startTime)
-
-    # print('Execution time in seconds: ' + str(executionTime))
-    # points = points.reshape(-1)
-    # write meshview json
-meshview.append({
-    "idx": str(1),
-    "count": str(len(temp_points)//3),
-    # "name"  : str(infoRow["name"].values[0]),
-    "name":"1",
-    "triplets": points_all,
-    # "r": str(infoRow['r'].values[0]),
-    # "g": str(infoRow['g'].values[0]),
-    # "b": str(infoRow['b'].values[0])
-    "r": 0,
-    "g": 0,
-    "b": 255
-})
-with open(f"colour_test_nonlin_notenpc.json", "w") as f:
-    json.dump(meshview, f)
\ No newline at end of file