Skip to content
Snippets Groups Projects
Commit 96bec45e authored by Andrei Mancu's avatar Andrei Mancu
Browse files

Update cycle detection

parent fe7742dc
No related branches found
No related tags found
No related merge requests found
......@@ -2,6 +2,12 @@ import networkx as nx
import numpy as np
import os
import time
import pandas as pd
from syconn import global_params
from syconn.reps.super_segmentation import SuperSegmentationDataset
from syconn.proc.meshes import calc_contact_syn_mesh, mesh2obj_file_colors, merge_meshes, gen_mesh_voxelmask
from scipy.spatial import cKDTree
from tqdm import tqdm
......@@ -9,17 +15,41 @@ from tqdm import tqdm
# from dgl.geometry import farthest_point_sampler
from networkx import from_edgelist, find_cycle, cycle_basis, simple_cycles, minimum_cycle_basis
from syconn import global_params
from syconn.reps.super_segmentation import SuperSegmentationDataset
from syconn.proc.meshes import calc_contact_syn_mesh, mesh2obj_file_colors, merge_meshes, gen_mesh_voxelmask
PINK = np.array([10., 255., 10., 255.])
BLUE = np.array([255., 125., 125., 255.])
GREY = np.array([180., 180., 180., 255.])
RED = np.array([10., 10., 255., 255.])
wd = '/cajal/nvmescratch/projects/from_ssdscratch/songbird/j0251/j0251_72_seg_20210127_agglo2/'
# wd = '/cajal/nvmescratch/projects/from_ssdscratch/songbird/j0251/j0251_72_seg_20210127_agglo2/'
wd = '/cajal/nvmescratch/projects/data/songbird/j0251/j0251_72_seg_20210127_agglo2_syn_20220811_celltypes_20230822/'
global_params.wd = wd
ids_dict = {}
file_name = '/cajal/nvmescratch/users/amancu/merge-error-detection/merge-error-detection/alex_self_merge_cycles_all.csv'
def recursive_find(edge_tuple, edge_list, cur_depth, min_depth, max_depth):
# Maximum depth was achieved
if cur_depth == max_depth:
return 1
# Prepare elements and pop form list
input_elem, output_elem = edge_tuple[0], edge_tuple[1]
cur_count = 0
# Run through all the branches of the input_node
for next_edge in edge_list:
if output_elem == next_edge[0]:
reach_max = recursive_find(next_edge, edge_list, cur_depth + 1, min_depth, max_depth)
cur_count += reach_max
if min_depth > cur_depth:
return cur_count
else:
if cur_count > 0:
return 1
else:
return 0
def find_longest_cycle(ssv):
print("skeleton_path:", ssv.skeleton_path)
......@@ -34,12 +64,13 @@ def find_longest_cycle(ssv):
scaled_nodes = np.array(scaled_nodes)
dists, idxs = tree.query(scaled_nodes, k=2)
dists, idxs = tree.query(scaled_nodes, k=2) # TODO : increase k
# Eliminate self nodes
idxs = np.array([elem[1] for elem in idxs])
dists = np.array([elem[1] for elem in dists])
candidate_idxs = np.where(dists < 150, True, False)
# Focus on close non-connected nodes in the graph
candidate_idxs = np.where(dists < 150, True, False)
add_edges = []
for src, target in enumerate(idxs):
......@@ -47,15 +78,19 @@ def find_longest_cycle(ssv):
continue
else:
path_length = nx.shortest_path_length(graph,src,target)
# Omit close non-connected nodes that are close to eachother in the skeleton graph (for ex: dendrite nodes)
# Add really close non-connected nodes that are far apart in the graph traversal to the add list
if path_length > 100:
add_edges.append([src,target])
# Distances for the src->target paths that correpsond to True in the candidate_idxs
graph.add_weighted_edges_from([(elem[0], elem[1], 1) for elem in add_edges])
# Generate cycles in the cell with the new graph nodes
cycles = cycle_basis(graph)
# Detect the largest cycle
if len(cycles) == 0:
return None
length = 0
......@@ -81,7 +116,46 @@ def find_longest_cycle(ssv):
cycle_candidates = longest_cycle[cycle_candidates_idx]
# TODO Exclude nodes inside soma
# TODO NOTE do we want to exclude candidate nodes inside soma??
# Pinpoint the merge spot by identifying nodes
# that have >3 leaf nodes in BFS with a 5<depth<15 nodes (to exclude dendrite paths)
depth_limit = 50 # BFS lookout limit
nr_branches = 4 # number of minimum branches expected throught the branching paths
min_lookout_depth = 5 # the max distance of the merge nodes such that a merge can occur
max_depth_dfs = 20 # maximum depth of the dfs for the 4 branches
min_depth_dfs = 7
# -> only clear axon/dendrite passings will be identified
# For example bouton/dendrite mergers will not be identified
merge_candidates = []
for merge_location_candidate in cycle_candidates:
bfs = nx.bfs_tree(graph, source = merge_location_candidate, depth_limit=depth_limit)
bfs_path_lengths = nx.bfs_layers(bfs, [merge_location_candidate])
lengths_list = np.array([len(lis) for lis in list(bfs_path_lengths)])
# All subsequent paths (>5 depth) should have at least 4 ramifications
res1 = np.all(lengths_list[min_lookout_depth:] >= nr_branches)
# dfs = nx.dfs_tree(graph, source = merge_location_candidate, depth_limit = max_depth_dfs)
dfs = nx.dfs_edges(graph, source=merge_location_candidate, depth_limit=max_depth_dfs)
list_dfs = list(dfs)
# Traverse each branch from the initial node
# all to all hops -> if two candidates are too close to eachother, regard them as individual branch
number_of_limit_reaches = 0
for next_edge in list_dfs:
if merge_location_candidate == next_edge[0]:
limit_reach = recursive_find(edge_tuple=next_edge, edge_list=list_dfs, cur_depth=1, min_depth=min_depth_dfs, max_depth=max_depth_dfs)
number_of_limit_reaches += limit_reach
res2 = False
if number_of_limit_reaches > 3:
res2 = True
if res1 and res2:
merge_candidates.append(merge_location_candidate)
if len(longest_cycle) != 0:
......@@ -93,7 +167,7 @@ def find_longest_cycle(ssv):
try:
np.put_along_axis(colors, mask, PINK, axis=0)
except:
print("No foreground labels in original context.")
print("No cycle nodes.")
pass
mask = np.array([[x] for x in cycle_candidates])
try:
......@@ -101,32 +175,84 @@ def find_longest_cycle(ssv):
except:
print("No foreground labels in original context.")
pass
mesh2obj_file_colors(os.path.expanduser(
f'test/cycle_{ssv.id}.ply'),
[np.array([]), scaled_nodes, np.array([])], colors)
mask = np.array([[x] for x in merge_candidates])
try:
np.put_along_axis(colors, mask, RED, axis=0)
except:
print("No foreground labels in original context.")
pass
mesh2obj_file_colors(os.path.expanduser(
f'alex_segmentation_all/cycle_{ssv.id}_loc_new.ply'),
[np.array([]), scaled_nodes / ssv.scaling, np.array([])], colors)
for merge_cand in merge_candidates:
print(f'merge cand: {merge_cand}')
ids_dict[f'{ssv.id}'] = [scaled_nodes[merge_cand] / ssv.scaling]
print(f"debug: found cycles, longest: {length}")
if length > 10:
# points = ssv.skeleton["nodes"][longest_cycle]
# three_idx = farthest_point_sampler(points, 3)
# three_pos = points[three_idx]
return length#, three_pos
return length
else:
return None
ssd = SuperSegmentationDataset(working_dir=global_params.config.working_dir)
ssv_ids = ssd.ssv_ids
cycles = []
# lis = [1587705582, 2525003196, 1453824753, 240638196, 50644658, 1251320138, 892693433, 1782961402, 1035245748, 1084271419, 4488813, 1581162187, 1253780676, 17986404, 620481014, 1981755000, 42752011, 1562064217, 1812246333, 1392116416, 2701401582, 1147566839, 38295053, 41639160, 1930845030, 583636574, 1720103625, 2000817469, 427017217, 1613922103, 1771722082, 47878065, 13159824, 1481931172, 2148059011, 30738782, 31786410, 73354653]
lis = [4726925, 824639693, 1299023959, 1126984089, 13327739, 331699998, 1957563793]
# lis = [824639693]
lis = sorted(lis)
# lis = [4726925, 824639693, 1299023959, 1126984089, 13327739, 331699998, 1957563793]
all_lis = [3374061, 3442140, 3948140, 4118195, 4150357, 4352555, 7560202, 7626258, 7962964, 8504301, 9885867, 11406245, 11572909, 12146425, 12585048, 15320718, 15420661, 16163954, 18962132, 19775310, 21897297, 23314382, 23788710, 24194587, 26114836, 26621052, 31278593, 38297008, 38297192, 39111453, 39782757, 39884461, 42111883, 43223315, 44874367, 44876588, 48855586, 49666859, 49769370, 50648088, 52229570, 52330826, 52332920, 52433603, 53378460, 55266237, 60330383, 61912928, 64376014, 67046128, 67551068, 74565413, 75476834, 75682236, 82192745, 82937782, 84555458, 86140676, 93431486, 126289248, 129934924, 163842826, 167456363, 196875835, 242660426, 266718802, 291719074, 333150399, 361963956, 390072183, 396750419, 432988156, 451680496, 481406285, 564472565, 662759023, 687119417, 756622970, 775989712, 821097737, 825519277, 883552299, 895023368, 980385678, 992703446, 1054782299, 1130665247, 1204864258, 1259113748, 1266095423, 1298013841, 1299023959, 1306618463, 1312459848, 1380947726, 1463980893, 1479772383, 1554336852, 1590608874, 1598028160, 1645435353, 1703808260, 1846729013, 1926928359, 1985032745, 2001761429, 2072584686, 2307245899, 4488613, 10157981, 24596473, 25372570, 26353030, 26958335, 28209513, 33605129, 45452841, 52904364, 78211969, 175077676, 244383709, 4150116, 5803292, 10190245, 17748051, 32019367, 32356701, 38767646, 40660209, 46092331, 50614045, 61710549, 122781754, 347218588, 354137851, 600576536, 606343972, 775215267, 824236472, 968947981, 1101410378, 1389045173, 1453824753, 1557170527, 1714433149, 1753707157, 1858603352, 1903488823, 1930165972, 2834364, 3341042, 4151286, 8266372, 11541121, 11675355, 13394835, 15388259, 16060305, 18963088, 20650554, 21022344, 23552465, 26252706, 26519659, 26790127, 27531858, 28814601, 33643243, 34178656, 39847003, 44671829, 46999900, 58336508, 58606391, 62419039, 63264864, 63336405, 63431281, 64210280, 67009484, 69200962, 72541034, 72879213, 83542452, 85598412, 207333487, 314561390, 328866567, 369961324, 379072583, 453097983, 479415063, 480798536, 580061165, 597605469, 686948959, 792864986, 809390056, 844007128, 946881951, 1075431186, 1179179749, 1193993147, 1431893661, 1594081077, 1724891551, 2094180358, 4352945, 4554974, 7928917, 15521116, 17106332, 20718752, 21697917, 22606703, 25409235, 30703527, 41095452, 42918146, 43865970, 44233162, 44941740, 48720742, 48926571, 51993379, 55772296, 62152255, 67750051, 88265889, 139382594, 173325651, 239792676, 398809969, 480462378, 481134355, 514467771, 639982340, 657459333, 756218606, 778925828, 1040205511, 1091356471, 1270685701, 1567090303, 1950679306, 3475290, 4084266, 4554921, 4555676, 4590645, 4726925, 11373887, 13733529, 15318237, 17415252, 22572226, 24530588, 29222156, 30298743, 35561999, 37251840, 38433620, 41972704, 49362703, 50376575, 51021023, 51353357, 54493399, 58306741, 59688447, 65830797, 69538197, 73452048, 75479167, 75546553, 76693645, 79998766, 177407557, 185201455, 232031927, 247147704, 283013933, 329642864, 448206654, 483970022, 504011358, 523544385, 566902791, 617645784, 642210125, 647000676, 702065851, 731216132, 731787776, 742619192, 747983832, 760129356, 804565781, 826024820, 843906326, 864117497, 870732015, 928023064, 945465952, 945905333, 1027589271, 1080627023, 1095370130, 1208638681, 1237215610, 1244030544, 1284722883, 1288095636, 1292077191, 1300075738, 1387828451, 1403112515, 1425450058, 1439417312, 1457569565, 1474000401, 1481427275, 1497316472, 1511049548, 1574448572, 1582074009, 1680660653, 1739266459, 1782961402, 1786233482, 1819634696, 1872809754, 1885902289, 2144784852, 2164997307, 2209030496, 2221815992, 2223572564, 2224684990, 2289634586, 2310381204, 2800426, 3104471, 3476114, 4219934, 4423089, 4457066, 7562545, 10360499, 12417788, 12720603, 13327739, 20886264, 25035193, 34856135, 35531229, 42042881, 43524762, 44271205, 46191576, 47910863, 48723488, 54256546, 56514562, 66000105, 76625102, 79930190, 81890441, 82967407, 85903087, 126798179, 135537053, 137592824, 243603620, 252241982, 260711535, 316179505, 350056908, 366787962, 379173238, 414160769, 441220870, 460321100, 481269892, 562314346, 624496225, 648856098, 653344899, 657494197, 691200729, 696261025, 743802111, 777978626, 821198622, 828725534, 839082899, 883348491, 884326921, 896472954, 924915907, 988549701, 993277627, 1034372879, 1084879573, 1107246316, 1155532413, 1213359231, 1234246033, 1305267601, 1313061357, 1362561086, 1370318073, 1395622088, 1458412290, 1481728632, 1489184340, 1495426731, 1586693935, 1605420115, 1721587518, 1723205375, 1736632955, 1741224471, 1765819411, 1842821366, 1843958990, 1906448471, 1944842013, 1951084693, 1962118669, 2006891760, 2017924639, 2052813803, 2137836140, 2177987581, 2207310366, 2207711447, 2261188431, 2365583723, 2382320636, 2426483413, 4117198, 4185991, 5469242, 10293994, 11709373, 12013841, 12382639, 15724767, 21224556, 21324985, 21931278, 24397945, 34179069, 35700843, 37215277, 40154408, 45991236, 51284828, 53342974, 54155520, 62151441, 62489778, 63532486, 65761121, 74870486, 85804145, 180884175, 230412443, 325492436, 339189930, 372794785, 468213671, 468515837, 490749885, 524051704, 635092005, 730237444, 777573524, 802272084, 858247487, 899643648, 1009402097, 1025968640, 1080458707, 1080557833, 1084372433, 1095911839, 1226315192, 1251758207, 1251924311, 1359488848, 1368735046, 1491714153, 1503628583, 1542359051, 1548500233, 1577617217, 1613044317, 1660482168, 1685085253, 1763525469, 1781781736, 2016473479, 2017855719, 2082266283, 2484617743]
# Creating the full list with the provided numbers
msn_lis = [1409019620, 821199296, 1412629286, 494125967, 28847776, 61915167, 1989649398, 2207913970, 63267019, 153350916,
471721862, 191610264, 39745958, 2147213777, 2077309462, 2329851856, 1763556396, 650914835, 1739163107, 2120529178,
514535349, 37822591, 116269788, 400495976, 6005889, 1384726225, 891345302, 50581003, 1309757255, 200451476,
1440160916, 1845173808, 25878817, 1147566839, 1030791524, 331464329, 1277639098, 1853914575, 681448747, 8674415,
1040576090, 6243837, 1488645700, 43659617, 531777405, 2224618017, 64376800, 2028485768, 48891830, 851260521,
2101867900, 1279657785, 1867779637, 1579003657, 873464831, 1587773682, 1908676891, 54763400, 32560586, 2034289573,
1267951828, 926536379, 28612788, 16163390, 2187938116, 60765851, 65290075, 1516178466, 579893062, 1149962101,
65052732, 1596107746, 2202988516, 125110677, 1765075591, 35395840, 1405915815, 782499239, 1999538773, 126323449,
1055221151, 1461010673, 2489239243, 939291245, 14879595, 2155212907, 69369529, 33942510, 63768694, 216511926,
49028772, 2411705007, 939997975, 1311276430, 782468175, 1547114210, 1383075004, 18995727, 56346823, 63228802,
2175224217, 736783860, 2391294989, 726727190, 56685128, 1612875677, 1928615319, 2166110290, 63163484, 26656085,
1109172466, 68188688, 1739603386, 237026761, 1188935439, 73688255, 689547469, 78715445, 83101631, 760029868,
2000686284, 1075026208, 2097115717, 1990530856, 217287679, 1702625081, 232368717, 2271480709, 793465442, 514400642,
1900576566, 675780312, 1357834589, 3407745, 34178689, 2594403947, 452893519, 490549501, 56548248, 1450182359,
721598843, 772006115, 2352423465, 1289783356, 45856538, 42752087, 13394836, 2612355690, 1744665042, 1379026140,
1784347845, 984230591, 867932422, 2235007897, 1160355315, 1253140118, 55839773, 11811035, 1869739137, 52669830,
2235917697, 1483013247, 15285578, 2105750450, 727877352, 75142197, 184559514, 656786944, 124368221, 1437358844,
49260459, 13834068, 2074978276, 54255086, 1660347083, 1398288836, 4388176, 2052709852, 283957119, 92955864,
251064152, 54223316, 127268495, 1846729206, 3577511, 1954559192, 45314997, 1653261011, 1060418590, 1030320715,
1900307086, 37586442, 1778269342, 28243772, 1441676701, 2359677820, 54663752, 242358283, 217793582, 2470208497,
42816104, 1410671107, 84248858, 1083124391, 1375886127, 2443183674, 3611292, 1991948634, 762999195, 1612708766,
131114911, 969727313, 999211951, 263547936, 74602209, 12585042, 2089219959, 2408297386, 2081087397, 839049320,
1431588268, 1956111165, 2307516196, 2029192883, 1444884566, 43187281, 1887722279, 17612882, 1780868791, 1741559892,
1362997356, 718966951, 1596375102, 2367606012, 987775203, 2122653861, 769646045, 928358934, 691334640, 9686641,
719473875, 1951085460, 2405565387, 7625363, 7457320, 2280894944, 401676736, 1690920382, 134522045, 1145306089,
1357835733, 957139005, 75479700, 1667772593, 1790618552, 575371134, 888175629, 9413522, 697981781, 246136479,
1525827215, 544497886, 4689868, 2310954733, 48825309, 1435033272, 548513296, 46801779, 1134578096, 108409967,
554958569, 27666854, 2000817469, 883552535, 11709047, 1359825618, 2232882013, 2383867855, 674702297, 1871930079,
1410671170, 1659269594, 687217692, 4488813, 535826054, 20480625, 1979192738, 491493101, 73182182, 1091760690,
791879706, 626316394, 2292199814, 1270379116, 125480144, 1207862649, 1442083152, 25440205, 30804770, 214826481,
48585630, 8570811, 1210730689, 1752931013, 1623401057, 27536753, 1094122533, 976607602, 68965449, 1211505322,
729394231, 26620933, 2075249394, 1511386273, 73890741, 1444310381, 1579006226, 81853342, 811347515, 23114881,
2062697082, 632323660, 2095493202, 1771623570, 2150184032, 2229675047, 14342416, 514502646, 64679677, 551889247,
33942472, 1908573325, 61476038, 1403081645, 1749689716, 694372316, 652061236, 935443282, 927580652, 438822964,
2070155090, 1827258980, 242760870, 1447820145, 6478149, 777642760, 1400551423, 2518694336, 8435027, 67648794,
58572939, 2363926688, 1488984047, 1069459614, 1346634119, 475366461, 1685149698, 1472415640, 970937346, 54523987,
1479673305, 1400179638, 220262653, 2383093398, 1700129316, 1907128291, 70314757, 67886005, 448408923, 21357545, 43659944, 11471875, 50241152, 35528313, 42984861, 932543277, 1810556186, 2330998697, 143667260, 2281975812, 3239603, 1357465274, 2055509469, 321478270, 26620952, 1888464470, 1889511922, 10057501, 903187271, 833414423, 66064223, 54051630, 1550050999, 1545091903, 2188820042, 1952671652, 1925275398, 35395330, 1169667909, 310408805, 11743771, 18221446, 2019950103, 138640773, 9483469, 1806105267, 647372573, 1565133047, 1928886488, 2093674606, 2210617606, 1268152262, 2358631096, 1269133286, 56075935, 148625778, 1654949046, 2094044782, 1382701471, 160537972, 384911484, 7696343, 80301216, 492573399, 1706809325, 66199873, 2249446561, 169312765, 1301082936, 1483989921, 1677592230, 1488983436, 766304738, 137933546, 154597836, 889083509, 20009567, 693427302, 4455891, 785738703, 36304357, 1233907772, 21593701, 350899113, 1702253709, 1640104117, 1340357151, 1459630950, 954001233, 1782555583, 2097788407, 1241971802, 2507662138, 2035098655, 1025767058, 2022781828, 846269031, 48788071, 91133390, 1249565418, 13598183, 1288804576, 1407567692, 135839398, 2262639799, 1539760783, 23148457, 701997441, 1967620056, 50373896, 664983081, 26249791, 2297931827, 13160357, 70147439, 1986210785, 47204006, 1054209794, 683575066, 10427714, 1445660360, 985210893, 72406155, 1365022426, 830682574, 937163742, 83203183, 1706198966, 483868931, 996140949, 198089964, 1690343276, 1376222854, 28783242, 531879202, 529956956, 2179877543, 1777662468, 895596819, 474723544]
lis = sorted(all_lis)
print(f'A total number of {len(lis)} will be computed')
for ssv_id in tqdm(lis): # tqdm(ssd.ssv_ids): todo: uncomment
cycle = find_longest_cycle(ssd.get_super_segmentation_object(ssv_id))
if cycle is not None:
cycles.append(cycle)
print("found cycle: ", cycle)
else:
print(f'Did not find cycle')
\ No newline at end of file
print(f'Did not find cycle')
df = pd.DataFrame(ids_dict)
df.to_csv(file_name)
print(f'A total of {len(cycles)} cells were found to have cycles within them.')
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment