Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
Filename: /home/eric/Code/project/mno/nutrimorph/2-morgane/brick/component/soma.py
Line # Mem usage Increment Occurences Line Contents
============================================================
300 189.8 MiB 189.8 MiB 1 def __HysterisisImage__(image: array_t, low: float, high: float) -> array_t:
301
302 # Find the image max value
303 189.8 MiB 0.0 MiB 1 max_image = image.max()
304 # Find the image minimum non zero value
305 230.2 MiB 40.5 MiB 1 nonzero_sites = image.nonzero()
306 236.8 MiB 6.5 MiB 1 nonzero_values = image[nonzero_sites]
307 237.0 MiB 0.3 MiB 1 min_image = nonzero_values.min()
308
309 # Adapt the hysteresis thresholds to the max and min of the image
310 237.0 MiB 0.0 MiB 1 low = low * (max_image - min_image) + min_image
311 237.0 MiB 0.0 MiB 1 high = high * (max_image - min_image) + min_image
312
313 # Perform hysteresis
314 248.9 MiB 11.9 MiB 1 result = fl_.apply_hysteresis_threshold(image, low, high)
315 248.9 MiB 0.0 MiB 1 result = result.astype(np_.int8, copy=False)
316
317 248.9 MiB 0.0 MiB 1 return result
Filename: /home/eric/Code/project/mno/nutrimorph/2-morgane/brick/component/soma.py
Line # Mem usage Increment Occurences Line Contents
============================================================
320 208.6 MiB 208.6 MiB 1 def __MorphologicalCleaning__(image: array_t, result: array_t, selem) -> array_t:
321
322 208.6 MiB 0.0 MiB 16 for dep in range(image.shape[0]):
323 208.6 MiB 0.0 MiB 15 result[dep, :, :] = mp_.closing(result[dep, :, :], selem)
324 208.6 MiB 0.0 MiB 15 result[dep, :, :] = mp_.opening(result[dep, :, :], selem)
325
326 208.6 MiB 0.0 MiB 1 return result
hist_step_curvature=None,
103 number_of_bins_curvature=None,
104 hist_bins_borders_curvature=None,
105 with_plot=None,
106 interactive=None,
107 in_parallel=None,
108 ) -> pd_.DataFrame():
109 """"""
110 164.4 MiB 0.0 MiB 1 soma_t = sm_.soma_t
111 164.4 MiB 0.0 MiB 1 extension_t = xt_.extension_t
112
113 164.4 MiB 0.0 MiB 1 print(f"STARTED: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}")
114 164.4 MiB 0.0 MiB 1 start_time = tm_.time()
115
116 # TODO add proper warning with warn module
117 # --- Images
118
119 # Find the name of the file, eq. to the biological tested condition number x
120 164.4 MiB 0.0 MiB 1 name_file = os_.path.basename(data_path)
121 164.4 MiB 0.0 MiB 1 print("FILE: ", name_file)
122 164.4 MiB 0.0 MiB 1 name_file = name_file.replace(".tif", "")
123
124 164.4 MiB 0.0 MiB 1 name_dir = os_.path.dirname(data_path)
125 164.4 MiB 0.0 MiB 1 print("DIR: ", name_dir)
126
127 # Find the dimension of the image voxels in micron
128 # TODO do something more general, here too specific to one microscope metadata
129 164.9 MiB 0.5 MiB 2 size_voxel_in_micron = brick.image.unit.FindVoxelDimensionInMicron(
130 164.4 MiB 0.0 MiB 1 data_path, size_voxel_in_micron=size_voxel_in_micron
131 )
132
133 # Read the image
134 181.3 MiB 16.4 MiB 1 image = io_.volread(data_path)
135
136 # Image size verification - simple version without user interface
137 181.3 MiB 0.0 MiB 1 image = in_.ImageVerification(image, channel)
138
139 # iv_.image_verification(image, channel) # -> PySide2 user interface # TODO: must return the modified image!
140 # /!\ conflicts between some versions of PySide2 and Python3
141
142 # image = image[:, 800:, 300:]
143 181.3 MiB 0.0 MiB 1 image = image[:, 800:, 500:]
144 # image = image[:, 600:, 500:]
145 181.3 MiB 0.0 MiB 1 img_shape = image.shape
146 #
147 181.3 MiB 0.0 MiB 2 print(
148 181.3 MiB 0.0 MiB 1 f"IMAGE: s.{img_shape} t.{image.dtype} m.{np_.amin(image)} M.{np_.amax(image)}"
149 )
150
151 # Intensity relative normalization (between 0 and 1).
152 # Image for soma normalized for hysteresis
153 188.2 MiB 6.9 MiB 1 image_for_soma = in_.IntensityNormalizedImage(image)
154 # Image for extensions not normalized for frangi enhancement
155 189.8 MiB 1.5 MiB 1 image_for_ext = image.copy()
156
157 189.8 MiB 0.0 MiB 2 print(
158 189.8 MiB 0.0 MiB 1 f"NRM-IMG: t.{image_for_soma.dtype} m.{np_.amin(image_for_soma):.2f} M.{np_.amax(image_for_soma):.2f}"
159 )
160
161 # --- Initialization
162 189.8 MiB 0.0 MiB 1 som_nfo = {}
163 189.8 MiB 0.0 MiB 1 ext_nfo = {}
164 189.8 MiB 0.0 MiB 1 axes = {}
165
166 #
167 # --- Somas
168 189.8 MiB 0.0 MiB 1 print("\n--- Soma Detection")
169
170 # Change the soma parameters from micron to pixel dimensions, using the previously determined voxel size
171 189.8 MiB 0.0 MiB 2 soma_min_area_c = brick.image.unit.ToPixel(
172 189.8 MiB 0.0 MiB 1 soma_min_area_c, size_voxel_in_micron, dimension=(0, 1, 2)
173 )
174 189.8 MiB 0.0 MiB 2 soma_max_area_c = brick.image.unit.ToPixel(
175 189.8 MiB 0.0 MiB 1 soma_max_area_c, size_voxel_in_micron, dimension=(0, 1, 2)
176 )
177 189.8 MiB 0.0 MiB 2 soma_selem_c = mp_.disk(
178 189.8 MiB 0.0 MiB 1 brick.image.unit.ToPixel(soma_selem_micron_c, size_voxel_in_micron)
179 )
180
181 # - Create the maps for enhancing and detecting the soma
182 # Hysteresis thresholding and closing/opening
183 189.8 MiB 0.0 MiB 1 print("Hysteresis thresholding: ", end="")
184 208.6 MiB 18.9 MiB 2 som_nfo["map_small"] = soma_t.Map(
185 189.8 MiB 0.0 MiB 1 image_for_soma, soma_low_c, soma_high_c, soma_selem_c
186 )
187 208.6 MiB 0.0 MiB 1 print("Done")
188 # Deletion of too small elements
189 208.6 MiB 0.0 MiB 1 print("Morphological cleaning: ", end="")
190 225.2 MiB 16.6 MiB 2 som_nfo["map_small"], som_lmp_small = soma_t.DeleteSmallSomas(
191 208.6 MiB 0.0 MiB 1 som_nfo["map_small"], soma_min_area_c
192 )
193 # Deletion of too big elements (probably aggregation of cells)
194 # Separated from the deletion of small elements to avoid extensions to be detected where too big somas were just deleted
195 225.2 MiB 0.0 MiB 2 som_nfo["map"], som_lmp = soma_t.DeleteBigSomas(
196 225.2 MiB 0.0 MiB 1 som_nfo["map_small"], som_lmp_small, soma_max_area_c
197 )
198 225.2 MiB 0.0 MiB 1 print("Done")
199
200 # Labelling of somas, find the number of somas
201 225.2 MiB 0.0 MiB 1 print("Somas labelling: ", end="")
202 238.7 MiB 13.5 MiB 1 som_nfo["lmp"], n_somas = ms_.label(som_nfo["map"], return_num=True)
203 238.7 MiB 0.0 MiB 1 print("Done")
204
205 # Use relabel instead of label to optimize the algorithm. /!\ But BUG.
206 # som_nfo["lmp"] = relabel_sequential(som_lmp)[0]
207 # n_somas = som_nfo["lmp"].max()
208
209 # Create the distance and influence maps used next for connexions between somas-extensions and extensions-extensions
210 238.7 MiB 0.0 MiB 1 print("Distance and Influence Map creation: ", end="")
211 265.9 MiB 27.1 MiB 2 som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(
212 238.7 MiB 0.0 MiB 1 som_nfo["lmp"]
213 )
214 265.9 MiB 0.0 MiB 1 print("Done")
215
216 265.9 MiB 0.0 MiB 1 return
217
218 # Create the tuple somas, containing all the objects 'soma'
219 print("Somas creation: ", end="")
220 somas = tuple(
221 soma_t().FromMap(som_nfo["lmp"], uid) for uid in range(1, n_somas + 1)
222 )
223 print("Done")
224
225 print(f" n = {n_somas}")
226
227 elapsed_time = tm_.gmtime(tm_.time() - start_time)
228 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
229
230 if with_plot:
231 fb_.PlotSomas(somas, som_nfo, axes)
232 elif interactive:
233 window = smvl.soma_validation_window_t(
234 image_for_soma,
235 som_nfo["lmp"],
236 mip_axis=0,
237 gray_offset=0,
238 colormap="viridis",
239 )
240 n_somas = window.LaunchValidation()
241 print(f"Validated Somas: {n_somas}")
242
243 #
244 # -- Extentions
245 print("\n--- Extension Detection")
246
247 # Set to zeros the pixels belonging to the somas.
248 # Take into account the aggregation of cells detected as one big soma to avoid extension creation there
249 del_soma = (som_nfo["map_small"] > 0).nonzero()
250 image_for_ext[del_soma] = 0
251
252 # Change the extensions parameters from micron to pixel dimensions
253 ext_min_area_c = brick.image.unit.ToPixel(
254 ext_min_area_c, size_voxel_in_micron, dimension=(0, 1, 2)
255 )
256
257 if ext_selem_micron_c == 0:
258 ext_selem_pixel_c = None
259 else:
260 ext_selem_pixel_c = mp_.disk(
261 brick.image.unit.ToPixel(ext_selem_micron_c, size_voxel_in_micron)
262 )
263
264 # - Perform frangi feature enhancement (via python or c - faster - implementation)
265 enhanced_ext, ext_scales = extension_t.EnhancedForDetection(
266 image_for_ext,
267 scale_range,
268 scale_step,
269 alpha,
270 beta,
271 frangi_c,
272 bright_on_dark,
273 method,
274 diff_mode=diff_mode,
275 in_parallel=in_parallel,
276 )
277
278 elapsed_time = tm_.gmtime(tm_.time() - start_time)
279 print(f"Elapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}\n")
280
281 # - Creation of the enhanced maps
282 # Enhanced the contrast in frangi output
283 print("Enhancing Contrast: ", end="")
284 if adapt_hist_equalization:
285 # necessary but not sufficient, histogram equalisation (adaptative)
286 enhanced_ext = ex_.equalize_adapthist(enhanced_ext)
287 else:
288 # necessary but not sufficient, histogram equalisation (global)
289 enhanced_ext = ex_.equalize_hist(enhanced_ext)
290
291 # Rescale the image by stretching the intensities between 0 and 255, necessary but not sufficient
292 enhanced_ext = ex_.rescale_intensity(enhanced_ext, out_range=(0, 255))
293 print("Done")
294 # Hysteresis thersholding
295 print("Hysteresis Thresholding: ", end="")
296 ext_nfo["coarse_map"] = extension_t.CoarseMap(
297 enhanced_ext, ext_low_c, ext_high_c, ext_selem_pixel_c
298 )
299 print("Done")
300 # Deletion of too small elements
301 print("Morphological cleaning: ", end="")
302 ext_nfo["coarse_map"], ext_lmp = extension_t.FilteredCoarseMap(
303 ext_nfo["coarse_map"], ext_min_area_c
304 )
305 print("Done")
306 # Creation of the extensions skeleton
307 print("Skeletonization and Thinning: ", end="")
308 ext_nfo["map"] = extension_t.FineMapFromCoarseMap(ext_nfo["coarse_map"])
309 print("Done")
310 # Deletion of extensions found within the somas
311 print("Map Labelling: ", end="")
312 ext_nfo["map"][som_nfo["map"] > 0] = 0
313 # Labelling of extensions and number of extensions determined
314 ext_nfo["lmp"], n_extensions = ms_.label(ext_nfo["map"], return_num=True)
315 print("Done")
316
317 # Use relabel instead of label to optimize the algorithm. BUT PROBLEM WITH THE NUMBER OF EXTENSIONS DETECTED !
318 # ext_nfo["lmp"] = relabel_sequential(ext_lmp)[0]
319 # n_extensions = ext_nfo["lmp"].max()
320
321 # Create the tuple extensions, containing all the objects 'extension'
322 print("Extensions creation: ", end="")
323 extensions = tuple(
324 extension_t().FromMap(ext_nfo["lmp"], ext_scales, uid)
325 for uid in range(1, n_extensions + 1)
326 )
327 print("Done")
328 print(f" n = {n_extensions}")
329
330 # Create global end point map for extensions
331 glob_ext_ep_map = xt_.EndPointMap(ext_nfo["lmp"] > 0)
332 all_end_points = list(
333 zip(*glob_ext_ep_map.nonzero())
334 ) # updated in ValidateConnexion
335
336 #
337 elapsed_time = tm_.gmtime(tm_.time() - start_time)
338 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
339
340 if with_plot:
341 fb_.PlotExtensions(extensions, ext_nfo, img_shape)
342
343 # -- Preparation for Connections
344 # Calculate the COSTS of each pixel on the image and DILATE the existing extensions to
345 # avoid tangency of connexions to extensions
346 dijkstra_costs = dk_.DijkstraCosts(image, som_nfo["map"], ext_nfo["map"])
347
348 if dilatation_erosion:
349 # Dilate all extensions
350 for extension in extensions:
351 cn_.Dilate(extension.sites, dijkstra_costs)
352
353 # -- Soma-Extention
354 print("\n--- Soma <-> Extension")
355
356 # Change the extensions parameters from micron to pixel dimensions
357 max_straight_sq_dist_c = brick.image.unit.ToPixel(
358 max_straight_sq_dist_c, size_voxel_in_micron
359 )
360 max_weighted_length_c = brick.image.unit.ToPixel(
361 max_weighted_length_c, size_voxel_in_micron
362 )
363
364 # Find the candidate extensions, with their endpoints, for connexion to the somas
365 candidate_conn_nfo = cn_.CandidateConnections(
366 somas,
367 som_nfo["influence_map"],
368 som_nfo["dist_to_closest"],
369 extensions,
370 max_straight_sq_dist_c,
371 )
372
373 # For each candidate, verify that it is valid based on the shortest path and the maximum allowed straight distance
374 for ep_idx, soma, extension, end_point in candidate_conn_nfo:
375 # Only try to connect if the extension is not already connected
376 if extension.is_unconnected:
377 print(f" Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end="")
378 if dilatation_erosion:
379 # Erode the end_point of the extension in the costs map
380 cn_.Erode(
381 (end_point,),
382 dijkstra_costs,
383 extension,
384 extensions,
385 image,
386 all_end_points=all_end_points,
387 )
388
389 # Use of Dijkstra shortest weighted path
390 path, length = cn_.ShortestPathFromToN(
391 image=image,
392 point=end_point,
393 extension=extension,
394 costs=dijkstra_costs,
395 candidate_points_fct=soma.ContourPointsCloseTo,
396 max_straight_sq_dist=max_straight_sq_dist_c,
397 erode_path=False,
398 )
399
400 # Keep the connexion only if inferior to the allowed max weighted distance
401 if path.__len__() > 0:
402 # length in pixel and max_weighted_length_c too
403 if length <= max_weighted_length_c:
404 # Validate and update all the fields + dilate again the whole extension
405 cn_.ValidateConnection(
406 soma,
407 extension,
408 end_point,
409 path,
410 dijkstra_costs,
411 all_ep=all_end_points,
412 )
413 if dilatation_erosion:
414 # Dilate extensions
415 cn_.Dilate(extension.sites, dijkstra_costs)
416 print(": Made")
417 else:
418 if dilatation_erosion:
419 cn_.Dilate(extension.sites, dijkstra_costs)
420 print("")
421 else:
422 if dilatation_erosion:
423 cn_.Dilate(extension.sites, dijkstra_costs)
424 print("")
425
426 brick.io.console.feedback.PrintConnectedExtensions(extensions)
427
428 #
429 elapsed_time = tm_.gmtime(tm_.time() - start_time)
430 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
431
432 if with_plot:
433 fb_.PlotSomasWithExtensions(somas, som_nfo, "all")
434
435 # -- Extention-Extention
436 print("\n--- Extension <-> Extension")
437
438 should_look_for_connections = True
439
440 # Update the soma maps with next Ext-Ext connexions
441 while should_look_for_connections:
442
443 if dilatation_erosion:
444 # Update the costs by dilating everything again plus the contour points of the soma
445 for soma in somas:
446 cn_.Dilate(soma.contour_points, dijkstra_costs)
447 # Below probably not necessary but security
448 for extension in soma.extensions:
449 cn_.Dilate(extension.sites, dijkstra_costs)
450
451 # Update the final somas + ext map by adding the new extensions and their connexion paths
452 som_nfo["soma_w_ext_lmp"] = soma_t.SomasLMap(somas)
453
454 # Update the influence map
455 som_nfo["dist_to_closest"], som_nfo["influence_map"] = soma_t.InfluenceMaps(
456 som_nfo["soma_w_ext_lmp"]
457 )
458
459 # Find the candidate extensions for connexion to the primary extensions
460 candidate_conn_nfo = cn_.CandidateConnections(
461 somas,
462 som_nfo["influence_map"],
463 som_nfo["dist_to_closest"],
464 extensions,
465 max_straight_sq_dist_c,
466 )
467
468 should_look_for_connections = False
469
470 # For each candidate, verify if valid based on the shortest path and the maximum allowed straight distance
471 for ep_idx, soma, extension, end_point in candidate_conn_nfo:
472 # Only try to connect if the extension is not already connected
473 if extension.is_unconnected:
474 print(
475 f" Soma.{soma.uid} <-?-> Ext.{extension.uid}({ep_idx})", end=""
476 )
477 if dilatation_erosion:
478 # Erode the end_point of the extension in the costs map
479 cn_.Erode(
480 (end_point,),
481 dijkstra_costs,
482 extension,
483 extensions,
484 image,
485 all_end_points=all_end_points,
486 )
487 # Use of Dijkstra shortest weighted path
488 path, length = cn_.ShortestPathFromToN(
489 image=image,
490 point=end_point,
491 extension=extension,
492 costs=dijkstra_costs,
493 candidate_points_fct=soma.ExtensionPointsCloseTo,
494 extensions=extensions,
495 all_end_points=all_end_points,
496 max_straight_sq_dist=max_straight_sq_dist_c,
497 erode_path=dilatation_erosion,
498 )
499
500 # Keep the connexion only if inferior to the allowed max weighted distance
501 if path.__len__() > 0:
502 if length <= max_weighted_length_c:
503 # Verify that the last element of the connexion path is contained in the extension to be connected
504 # If so, returns the extensions containing the path last site.
505 tgt_extenstion = extension_t.ExtensionContainingSite(
506 extensions, path[-1]
507 )
508 # Validate connexion: update the glial component objects and store more info in their fields
509 cn_.ValidateConnection(
510 tgt_extenstion,
511 extension,
512 end_point,
513 path,
514 dijkstra_costs,
515 all_ep=all_end_points,
516 )
517 if dilatation_erosion:
518 # Dilate extensions
519 cn_.Dilate(extension.sites, dijkstra_costs)
520
521 should_look_for_connections = True
522
523 print(": Made")
524 else:
525 if dilatation_erosion:
526 cn_.Dilate(extension.sites, dijkstra_costs)
527 print("")
528 else:
529 if dilatation_erosion:
530 cn_.Dilate(extension.sites, dijkstra_costs)
531 print("")
532
533 # -- Create an extension map containing all the ext + connexions, without the somas.
534 # Used for the graphs extraction
535 ext_nfo["lmp_soma"] = som_nfo["soma_w_ext_lmp"].copy()
536 ext_nfo["lmp_soma"][som_nfo["map"] > 0] = 0
537 #
538
539 brick.io.console.feedback.PrintConnectedExtensions(extensions)
540
541 #
542 elapsed_time = tm_.gmtime(tm_.time() - start_time)
543 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
544
545 if with_plot:
546 fb_.PlotSomasWithExtensions(somas, som_nfo, "with_ext_of_ext")
547
548 # -- Summary
549 print("\n")
550 for soma in somas:
551 print(soma)
552
553 elapsed_time = tm_.gmtime(tm_.time() - start_time)
554 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
555
556 if with_plot:
557 pl_.show()
558
559 if save_images is not None:
560 po_.MaximumIntensityProjectionZ(
561 som_nfo["soma_w_ext_lmp"],
562 block=False,
563 output_image_file_name=str(pathlib.Path(save_images) / f"MIP_{name_file}.png"),
564 )
565
566 # --- Extract all the extensions of all somas as a graph
567 print("\n--- Graph extraction")
568
569 # Create the graphs with SKLGraph module (available on Eric Debreuve Gitlab)
570 for soma in somas:
571 print(f" Soma {soma.uid}", end="")
572 # Create SKLGraph skeletonized map
573 ext_map = skl_map_t.FromShapeMap(
574 ext_nfo["lmp_soma"] == soma.uid,
575 store_widths=True,
576 skeletonize=True,
577 do_post_thinning=True,
578 )
579
580 # Create the graph from the SKLGaph skeletonized map
581 soma.skl_graph = skl_graph_t.FromSkeleton(
582 ext_map, end_point, size_voxel=size_voxel_in_micron
583 )
584
585 # --- Find the root of the {ext+conn} graphs.
586 # Roots are the nodes of degree 1 that are to be linked to the soma
587 soma.graph_roots = ge_.FindGraphsRootWithNodes(soma)
588
589 # Add a node "soma" and link it to the root nodes
590 soma_node = (
591 f"S-{int(soma.centroid[0])}-{int(soma.centroid[1])}-{int(soma.centroid[2])}"
592 )
593 soma.skl_graph.add_node(soma_node, soma=True, soma_nfo=soma)
594
595 for node in soma.graph_roots.values():
596 soma.skl_graph.add_edge(node, soma_node, root=True)
597
598 if save_images is not None:
599 nx_.draw_networkx(soma.skl_graph)
600 pl_.savefig(str(pathlib.Path(save_images) / f"graph_{name_file}_soma{soma.uid}.png"))
601 # pl_.show(block=True)
602 pl_.close()
603
604 print(": Done")
605
606 elapsed_time = tm_.gmtime(tm_.time() - start_time)
607 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
608
609 # --- Extract features
610 print("\n--- Features Extraction\n")
611
612 # Parameters
613 if hist_bins_borders_length is None:
614 number_of_bins_length = int(number_of_bins_length)
615 bins_length = np_.linspace(
616 hist_min_length,
617 hist_min_length + hist_step_length * number_of_bins_length,
618 num=number_of_bins_length,
619 )
620 bins_length[-1] = np_.inf
621 else:
622 bins_length = np_.array(hist_bins_borders_length)
623 bins_length[-1] = np_.inf
624
625 if hist_bins_borders_curvature is None:
626 number_of_bins_curvature = int(number_of_bins_curvature)
627 bins_curvature = np_.linspace(
628 hist_min_curvature,
629 hist_min_curvature + hist_step_curvature * number_of_bins_curvature,
630 num=number_of_bins_curvature,
631 )
632 bins_curvature[-1] = np_.inf
633 else:
634 bins_curvature = np_.array(hist_bins_borders_curvature)
635 bins_curvature[-1] = np_.inf
636
637 # Pandas dataframe creation with all the measured features
638 features_df = ge_.ExtractFeaturesInDF(
639 name_file,
640 somas,
641 size_voxel_in_micron,
642 bins_length,
643 bins_curvature,
644 ext_scales,
645 )
646
647 # Save the pandas df into .csv
648 if save_csv is not None:
649 features_df.to_csv(f"{save_csv}/{name_file}.csv")
650 else:
651 features_df.to_csv(f"{name_dir}/{name_file}.csv")
652
653 elapsed_time = tm_.gmtime(tm_.time() - start_time)
654 print(f"\nElapsed Time={tm_.strftime('%Hh %Mm %Ss', elapsed_time)}")
655 print(f"DONE: {tm_.strftime('%a, %b %d %Y @ %H:%M:%S')}\n")
656
657 return features_df
Filename: /home/eric/Code/project/mno/nutrimorph/2-morgane/nutrimorph.py
Line # Mem usage Increment Occurences Line Contents
============================================================
660 164.4 MiB 164.4 MiB 1 def Main():
661 """"""
662 164.4 MiB 0.0 MiB 1 if sy_.argv.__len__() < 2:
663 print("Missing parameter file argument")
664 sy_.exit(0)
665 164.4 MiB 0.0 MiB 1 argument = sy_.argv[1]
666 164.4 MiB 0.0 MiB 3 if not (
667 164.4 MiB 0.0 MiB 1 os_.path.isfile(argument)
668 164.4 MiB 0.0 MiB 1 and os_.access(argument, os_.R_OK)
669 164.4 MiB 0.0 MiB 1 and (os_.path.splitext(argument)[1].lower() == ".ini")
670 ):
671 print("Wrong parameter file path or type, or parameter file unreadable")
672 sy_.exit(0)
673
674 164.4 MiB 0.0 MiB 1 (
675 164.4 MiB 0.0 MiB 1 data_path,
676 164.4 MiB 0.0 MiB 1 save_images,
677 164.4 MiB 0.0 MiB 1 save_csv,
678 164.4 MiB 0.0 MiB 1 save_features_analysis,
679 164.4 MiB 0.0 MiB 1 statistical_analysis,
680 164.4 MiB 0.0 MiB 1 arguments,
681 164.4 MiB 0.0 MiB 1 ) = NutrimorphArguments(argument)
682
683 # Differentiate between path to a tiff file or to a repository
684 164.4 MiB 0.0 MiB 1 if pathlib.Path(data_path).is_file():
685 # Perform NutriMorph algorithm on the file entered in parameters
686 164.4 MiB 0.0 MiB 2 print(
687 164.4 MiB 0.0 MiB 1 "WARNING: Will not perform features analysis on a single image.\n For features analysis, "
688 "give a directory path.\n"
689 )
690 197.2 MiB 32.8 MiB 1 _ = NutriMorph(data_path, save_images, save_csv, **arguments)
691
692 elif pathlib.Path(data_path).is_dir():
693 # Keep the directory to the repository
694 # Initialize the future concatenated features
695 concatenated_features_df = pd_.DataFrame()
696 # Find all the tiff files in the parent and child repositories
697 for path in pathlib.Path(data_path).glob("**/*.tif"):
698 if path.is_file():
699 # Perform NutriMorph algorithm
700 features_df = NutriMorph(str(path), save_images, save_csv, **arguments)
701 # Concatenate all the dataframes
702 concatenated_features_df = concatenated_features_df.append(features_df)
703 # If some rows (ie. somas0) have NaN, drop them
704 # -- due to best fitting ellipsoid algo (JTJ is singular due to convex hull being flat)
705 concatenated_features_df = concatenated_features_df.dropna(
706 axis=0, how="any"
707 )
708
709 # Save to .csv in the parent repository
710 if save_csv is None:
711 concatenated_features_df.to_csv(f"{data_path}/features.csv")
712 else:
713 concatenated_features_df.to_csv(f"{save_csv}/features.csv")
714
715 # Clustering with this df and module features_analysis.py
716 if statistical_analysis:
717 os_.system(
718 f"feature_analysis.py {save_csv}/features.csv {save_features_analysis}"
719 )
720 else:
721 raise ImportError(f"{data_path}: Not a valid data path!")