1# -*- coding: utf-8 -*-
2"""
3Created on Wed Apr 5 09:45:56 2023
4
5@author: Bart Steemans. Govers Lab.
6"""
7from . import utilities as u
8import numpy as np
9import logging
10
11bactoscoop_logger = logging.getLogger("logger")
12
13
14class Features:
15
16 def __init__(self, image_obj):
17
18 self.image_obj = image_obj
19
20 def discard_large_meshes(self, cell, max_mesh_size=600):
21 """
22 Discard cells with meshes that are too small or large.
23
24 This method checks the size of a cell's mesh and flags it for discard if the mesh size does not fall within the acceptable range.
25
26 Parameters
27 ----------
28 cell : Cell
29 The Cell object representing the cell to be checked.
30
31 Returns
32 -------
33 None
34 """
35
36 if not 4 <= len(cell.mesh) <= max_mesh_size:
37
38 bactoscoop_logger.debug(
39 f"The mesh of Cell number {cell.cell_id} is too small or large"
40 )
41 cell.discard = True
42 def discard_bad_meshes(self, cell):
43
44 px = self.image_obj.px
45 step_length = u.get_step_length(cell.x1, cell.y1, cell.x2, cell.y2, px)
46 _, width_not_ordered = u.get_avg_width(cell.x1, cell.y1, cell.x2, cell.y2, px)
47 cell_area = u.get_area(cell.contour, px)
48 mesh_area = u.get_mesh_area(width_not_ordered, step_length)
49 ratio = mesh_area / cell_area
50
51 if 0.9 <= ratio <= 1.1:
52 cell.discard = False
53 else:
54 bactoscoop_logger.debug(f"Cell mesh with cell id {cell.cell_id} is discarded due to bad mesh quality")
55 cell.discard = True
56
57 def contour_selection(self, cell, use_shifted_contours):
58 if use_shifted_contours and cell.shifted_contour is not None:
59 return cell.shifted_contour, cell.shifted_midline, cell.shifted_mesh
60 else:
61 return cell.contour, cell.midline, cell.mesh
62
63 def morphological(self, cell, channel, all_data, *args):
64 """
65 Calculate morphological features for a cell.
66
67 """
68 cell.morphological_features = {}
69
70 # Prepare cropped mask and image for morphological measurements and constriction degree
71 if self.image_obj.inverted_image is None:
72 self.image_obj.get_inverted_image()
73 cell_mask, cropped_signal, _, _, x_offset, y_offset = u.get_cell_mask(
74 cell.contour, cell.midline, self.image_obj.inverted_image
75 )
76 cropped_contour, cropped_ml, cropped_msh = u.get_cropped_cell_data(
77 cell.contour, cell.midline, cell.mesh, x_offset, y_offset
78 )
79 cropped_signal_interp2d = u.interp2d(cropped_signal)
80
81 # Morphological measurements
82 px = self.image_obj.px
83 step_length = u.get_step_length(cell.x1, cell.y1, cell.x2, cell.y2, px)
84
85 cell.morphological_features["cell_length"] = u.get_length(step_length)
86 cell.morphological_features["cell_width"], width_not_ordered = u.get_avg_width(
87 cell.x1, cell.y1, cell.x2, cell.y2, px
88 )
89 cell.morphological_features["cell_aspect_ratio"] = (
90 cell.morphological_features["cell_length"]
91 / cell.morphological_features["cell_width"]
92 )
93 cell.profiling_data["cell_widthno"] = width_not_ordered
94
95 cell_area = u.get_area(cell.contour, px)
96 cell.morphological_features["cell_area"] = cell_area
97 cell.morphological_features["cell_volume"] = u.get_volume(
98 cell.x1, cell.y1, cell.x2, cell.y2, step_length, px
99 )
100 cell.morphological_features["cell_surface_area"] = u.get_surface_area(
101 width_not_ordered, step_length
102 )
103 cell.morphological_features["cell_SOV"] = u.get_surface_area_over_volume(
104 cell.morphological_features["cell_surface_area"],
105 cell.morphological_features["cell_volume"],
106 )
107
108 per, circ, poa, sph = u.get_cell_perimeter_measurements(
109 cell.contour, cell.morphological_features["cell_area"], px
110 )
111 cell.morphological_features["cell_perimeter"] = per
112 cell.morphological_features["cell_circularity"] = circ
113 cell.morphological_features["cell_compactness"] = poa
114 cell.morphological_features["cell_sphericity"] = sph
115 cell.morphological_features["cell_width_variability"] = (
116 u.get_cell_width_variability(width_not_ordered)
117 )
118
119 ecc, sol, conv = u.get_cell_convexity(cell_mask)
120 cell.morphological_features["cell_eccentricity"] = ecc
121 cell.morphological_features["cell_solidity"] = sol
122 cell.morphological_features["cell_convexity"] = conv
123
124 con_curvature, con_max_c, con_min_c, con_mean_c, con_std_c = (
125 u.get_curvature_characteristics(cell.contour, px)
126 )
127 cell.morphological_features["cell_max_contour_curvature"] = con_max_c
128 cell.morphological_features["cell_min_contour_curvature"] = con_min_c
129 cell.morphological_features["cell_mean_contour_curvature"] = con_mean_c
130 cell.morphological_features["cell_std_contour_curvature"] = con_std_c
131 cell.morphological_features["cell_contour_bending_energy"] = u.bending_energy(
132 cell.contour, con_curvature, px
133 )
134
135 mid_curvature, mid_max_c, mid_min_c, mid_mean_c, mid_std_c = (
136 u.get_curvature_characteristics(cell.midline, px)
137 )
138 cell.morphological_features["cell_max_midline_curvature"] = mid_max_c
139 cell.morphological_features["cell_min_midline_curvature"] = mid_min_c
140 cell.morphological_features["cell_mean_midline_curvature"] = mid_mean_c
141 cell.morphological_features["cell_std_midline_curvature"] = mid_std_c
142 cell.morphological_features["cell_midline_bending_energy"] = u.bending_energy(
143 cell.midline, mid_curvature, px
144 )
145 cell.morphological_features["cell_midline_sinuosity"] = u.sinuosity(
146 cell.midline
147 )
148
149 cell.morphological_features["cell_area_asymmetry"] = u.get_area_asymmetry(
150 step_length, width_not_ordered
151 )
152
153 regionprops_features = u.get_additional_regionprops_features(cell_mask, px)
154
155 # Store additional regionprops features
156 cell.morphological_features["max_feret_diameter"] = regionprops_features[
157 "F_MAX_FERET_DIAMETER"
158 ]
159 cell.morphological_features["equivalent_diameter"] = regionprops_features[
160 "F_EQUIVALENT_DIAMETER"
161 ]
162 cell.morphological_features["max_radius"] = regionprops_features[
163 "F_MAXIMUM_RADIUS"
164 ]
165 cell.morphological_features["mean_radius"] = regionprops_features[
166 "F_MEAN_RADIUS"
167 ]
168 cell.morphological_features["median_radius"] = regionprops_features[
169 "F_MEDIAN_RADIUS"
170 ]
171 cell.morphological_features["extent"] = regionprops_features["F_EXTENT"]
172
173 # Store moment-related features
174 moment_types = [
175 "moments",
176 "moments_central",
177 "moments_normalized",
178 "moments_hu",
179 "inertia_tensor",
180 "inertia_tensor_eigvals",
181 ]
182 for moment_type in moment_types:
183 for key in regionprops_features:
184 if f"{moment_type}" in key:
185 cell.morphological_features[key] = regionprops_features[key]
186
187 # Constriction degree calculations
188 width_px = u.get_width(cell.x1, cell.y1, cell.x2, cell.y2) * px
189
190 width_profile = width_px
191 (
192 cell.morphological_features["cell_constriction_degree_width"],
193 cell.morphological_features["cell_rel_constr_pos_width"],
194 constrDegree_abs,
195 ctpos,
196 _,
197 ) = u.constr_degree_single_cell_min(width_profile, step_length)
198
199 axial_profile = u.measure_smoothened_intensity(
200 cropped_ml, cropped_signal_interp2d, width=5
201 )
202 (
203 cell.morphological_features["cell_constriction_degree_axial"],
204 cell.morphological_features["cell_rel_constr_pos_axial"],
205 _,
206 _,
207 _,
208 ) = u.constr_degree_single_cell_min(axial_profile, step_length)
209
210 weighted_axial_profile = u.get_weighted_intprofile(axial_profile, width_px)
211 (
212 cell.morphological_features["cell_constriction_degree_width_axial"],
213 cell.morphological_features["cell_rel_constr_pos_width_axial"],
214 _,
215 _,
216 _,
217 ) = u.constr_degree_single_cell_min(weighted_axial_profile, step_length)
218
219 try:
220 cell.profile_mesh = u.get_profile_mesh(
221 cropped_msh, u.get_width(cell.x1, cell.y1, cell.x2, cell.y2)
222 )
223 except MemoryError:
224 cell.discard = True
225 pass
226 mesh_intensity = cropped_signal_interp2d.ev(
227 cell.profile_mesh[0], cell.profile_mesh[1]
228 ).T
229 mesh_intensity_profile = u.gaussian_smoothing(mesh_intensity)
230 (
231 cell.morphological_features["cell_constriction_degree_mesh"],
232 cell.morphological_features["cell_rel_constr_pos_mesh"],
233 _,
234 _,
235 _,
236 ) = u.constr_degree_single_cell_min(mesh_intensity_profile, step_length)
237
238 # Discarding cells with an erroneous mesh
239 self.discard_bad_meshes(cell)
240
241 def profiling(
242 self, cell, channel, all_data, use_shifted_contours, shift_signal, *args
243 ):
244 """
245 Calculate profiling features for a cell.
246
247 """
248 cell.profiling_features = {}
249 px = self.image_obj.px
250
251 if self.image_obj.channels is None or channel not in self.image_obj.channels:
252 raise AttributeError(
253 f"No objects are detected for channel {channel} or the channel is not loaded"
254 )
255 return
256
257 if channel not in self.image_obj.bg_channels:
258 self.image_obj.get_bg_channel(channel)
259 bactoscoop_logger.debug("Background subtraction completed")
260
261 contour, midline, mesh = self.contour_selection(cell, use_shifted_contours)
262
263 if shift_signal:
264 cropped_signal, cropped_mask, _, x_offset, y_offset = u.crop_image(
265 self.image_obj.bg_channels[channel],
266 contour,
267 phase_img=self.image_obj.image,
268 )
269 else:
270 cropped_signal, cropped_mask, _, x_offset, y_offset = u.crop_image(
271 self.image_obj.bg_channels[channel], contour, phase_img=None
272 )
273
274 cropped_contour, cropped_ml, cropped_msh = u.get_cropped_cell_data(
275 contour, midline, mesh, x_offset, y_offset
276 )
277 cropped_signal_interp2d = u.interp2d(cropped_signal)
278
279 try:
280 profile_mesh = u.get_profile_mesh(
281 cropped_msh, u.get_width(cell.x1, cell.y1, cell.x2, cell.y2)
282 )
283 except MemoryError:
284 cell.discard = True
285 pass
286
287 # Intensity features
288 cell.profiling_features["axial_intensity"] = u.measure_smoothened_intensity(
289 cropped_ml, cropped_signal_interp2d, width=10
290 )
291 cell.profiling_features["step_length_demograph"] = (
292 u.get_cumulative_step_length_midline(cropped_ml, px)
293 )
294 cell.profiling_features["normalized_axial_intensity"] = u.normalize_per_cell(
295 cell.profiling_features["axial_intensity"]
296 )
297 cell.profiling_features["cell_length"] = u.line_length(cropped_ml) * px
298 mesh_intensity = cropped_signal_interp2d.ev(profile_mesh[0], profile_mesh[1]).T
299 cell.profiling_features["raw_average_mesh_intensity"] = np.average(
300 mesh_intensity, axis=0
301 )
302 cell.profiling_features["average_mesh_intensity"] = u.gaussian_smoothing(
303 mesh_intensity
304 )
305 cell.profiling_features["raw_normalized_average_mesh_intensity"] = (
306 u.normalize_per_cell(cell.profiling_features["raw_average_mesh_intensity"])
307 )
308 cell.profiling_features["normalized_average_mesh_intensity"] = (
309 u.normalize_per_cell(cell.profiling_features["average_mesh_intensity"])
310 )
311
312 cell.profiling_features["normalized_average_mesh_intensity"] = (
313 u.normalize_per_cell(cell.profiling_features["average_mesh_intensity"])
314 )
315
316 step_length = u.get_step_length(cell.x1, cell.y1, cell.x2, cell.y2, px)
317 # useful for DAPI signal
318 (
319 cell.profiling_features["signal_constriction_degree_mesh"],
320 cell.profiling_features["signal_rel_constr_pos_mesh"],
321 _,
322 _,
323 _,
324 ) = u.constr_degree_single_cell_min(
325 cell.profiling_features["average_mesh_intensity"], step_length
326 )
327
328 (
329 glcm_dissimilarity,
330 glcm_correlation,
331 glcm_homogeneity,
332 glcm_energy,
333 glcm_contrast,
334 glcm_entropy,
335 ) = u.calculate_haralick_features(
336 cropped_contour,
337 cropped_signal,
338 )
339 cell.profiling_features["cell_mesh_intensity_mean"] = np.mean(
340 cell.profiling_features["average_mesh_intensity"]
341 )
342 cell.profiling_features["cell_mesh_intensity_std"] = np.std(
343 cell.profiling_features["average_mesh_intensity"]
344 )
345 cell.profiling_features["cell_mesh_kurtosis"] = u.get_kurtosis(
346 cell.profiling_features["average_mesh_intensity"]
347 )
348 cell.profiling_features["cell_mesh_skewness"] = u.get_skew(
349 cell.profiling_features["average_mesh_intensity"]
350 )
351
352
353
354 cell.profiling_features["cell_glcm_dissimilarity"] = glcm_dissimilarity
355 cell.profiling_features["cell_glcm_correlation"] = glcm_correlation
356 cell.profiling_features["cell_glcm_homogeneity"] = glcm_homogeneity
357 cell.profiling_features["cell_glcm_energy"] = glcm_energy
358 cell.profiling_features["cell_glcm_contrast"] = glcm_contrast
359 cell.profiling_features["cell_glcm_shannon_entropy"] = glcm_entropy
360
361 avg_width, _ = u.get_avg_width_no_px(cell.x1, cell.y1, cell.x2, cell.y2)
362 cell.profiling_features["radial_intensity_distribution"] = u.radial_intensities(
363 cropped_signal_interp2d,
364 cropped_contour,
365 width=avg_width,
366 min_val=np.min(cropped_signal),
367 max_val=np.max(cropped_signal),
368 num_erosions=10,
369 erosion_scale=None,
370 )
371
372 # Discarding cells with an erroneous mesh
373 self.discard_bad_meshes(cell)
374
375 def svm(self, cell, channel, all_data, *args):
376 """
377 Calculate features of a cell to be used in training the Support Vector Machine (SVM) or curate a dataset using a trained SVM model.
378
379 """
380 cell.svm_features = {}
381
382 # Preparation cropped mask and images for intensity measurements
383 if self.image_obj.inverted_image is None:
384 self.image_obj.get_inverted_image()
385 cropped_signal, cropped_mask, cropped_contour, x_offset, y_offset = (
386 u.crop_image(self.image_obj.inverted_image, cell.contour)
387 )
388 cropped_signal_interp2d = u.interp2d(cropped_signal)
389
390 px = self.image_obj.px
391 step_length = u.get_step_length(cell.x1, cell.y1, cell.x2, cell.y2, px)
392
393 # Contour is added here for labeling purposes
394 cell.svm_features["contour"] = cell.contour
395
396 # Morphological features
397 cell.svm_features["cell_length"] = u.get_length(step_length)
398 cell.svm_features["cell_width"], width_not_ordered = u.get_avg_width(
399 cell.x1, cell.y1, cell.x2, cell.y2, px
400 )
401 cell.svm_features["cell_area"] = u.get_area(cell.contour, px)
402 cell.svm_features["cell_volume"] = u.get_volume(
403 cell.x1, cell.y1, cell.x2, cell.y2, step_length, px
404 )
405 cell.svm_features["cell_surface_area"] = u.get_surface_area(
406 width_not_ordered, step_length
407 )
408 (
409 curvature,
410 cell.svm_features["max_curvature"],
411 cell.svm_features["min_curvature"],
412 cell.svm_features["mean_curvature"],
413 std_c,
414 ) = u.get_curvature_characteristics(cell.contour, px)
415 (
416 cell.svm_features["cell_perimeter"],
417 cell.svm_features["cell_circularity"],
418 cell.svm_features["cell_POA"],
419 cell.svm_features["cell_sphericity"],
420 ) = u.get_cell_perimeter_measurements(
421 cell.contour, cell.svm_features["cell_area"], px
422 )
423
424 # Cell surface intensity features
425 total_int, max_int, mean_int = u.get_total_phaco_intensity(
426 cropped_contour, cropped_signal.shape, cropped_signal_interp2d
427 )
428 cell.svm_features["phaco_total_intensity"] = total_int
429 cell.svm_features["phaco_max_intensity"] = max_int
430 cell.svm_features["phaco_mean_intensity"] = mean_int
431
432 # Contour intensity features
433 phaco_contour_intensity = u.get_contour_intensity(
434 cropped_contour, cropped_signal_interp2d
435 )
436 phaco_contour_intensity_variability = u.measure_contour_variability(
437 phaco_contour_intensity
438 )
439
440 sorted_phaco_contour_intensity = sorted(phaco_contour_intensity)
441 cell.svm_features["phaco_contour_peaks"] = u.find_signal_peaks(
442 phaco_contour_intensity, max_int
443 )
444 cell.svm_features["phaco_max_contour_intensity"] = np.max(
445 phaco_contour_intensity
446 )
447 cell.svm_features["phaco_mean_contour_intensity"] = np.mean(
448 phaco_contour_intensity
449 )
450 cell.svm_features["phaco_min_contour_intensity"] = np.mean(
451 sorted_phaco_contour_intensity[:10]
452 )
453 cell.svm_features["phaco_max_contour_variability"] = np.max(
454 phaco_contour_intensity_variability
455 )
456 cell.svm_features["phaco_mean_contour_variability"] = np.mean(
457 phaco_contour_intensity_variability
458 )
459
460 # Midline intensity features
461 # If no midline available calculate it with mesh2midline( ).
462 cropped_midline = cell.midline.copy()
463 cropped_midline -= [x_offset, y_offset]
464 cell.profiling_data["phaco_axial_intensity"] = u.measure_smoothened_intensity(
465 cropped_midline, cropped_signal_interp2d, width=5
466 )
467
468 cell.svm_features["midline_kurtosis"] = u.get_kurtosis(
469 cell.profiling_data["phaco_axial_intensity"]
470 )
471 cell.svm_features["midline_skewness"] = u.get_skew(
472 cell.profiling_data["phaco_axial_intensity"]
473 )
474
475 # Expanded contour intensity features
476 expanded_contour = u.expand_contour(cropped_contour, scale=2)
477 eroded_contour = u.erode_contour(cropped_contour, scale=2)
478 cell.profiling_data["phaco_expanded_contour_intensity"] = (
479 u.get_contour_intensity(expanded_contour, cropped_signal_interp2d)
480 )
481 cell.profiling_data["phaco_eroded_contour_intensity"] = u.get_contour_intensity(
482 eroded_contour, cropped_signal_interp2d
483 )
484
485 cell.svm_features["phaco_max_expanded_contour_intensity"] = np.max(
486 cell.profiling_data["phaco_expanded_contour_intensity"]
487 )
488 cell.svm_features["phaco_mean_expanded_contour_intensity"] = np.mean(
489 cell.profiling_data["phaco_expanded_contour_intensity"]
490 )
491
492 # Mesh gradient features
493 cell.svm_features["phaco_cell_edge_gradient"] = np.average(
494 cell.profiling_data["phaco_eroded_contour_intensity"]
495 - cell.profiling_data["phaco_expanded_contour_intensity"]
496 )
497
498 # Discarding cells with an erroneous mesh
499 self.discard_bad_meshes(cell)
500
501 def objects(
502 self, cell, channel, all_data, use_shifted_contours, *args,
503 retain_contour_on_object_mesh_failure=False,
504 ):
505 """
506 Calculates features related to the identified objects within a cell.
507
508 """
509 cell.objects_features = {}
510
511 px = self.image_obj.px
512
513 if self.image_obj.channels is None or channel not in self.image_obj.channels:
514 raise AttributeError(
515 f"No objects are detected for channel {channel} or the channel is not loaded"
516 )
517 return
518
519 if channel not in self.image_obj.bg_channels:
520 self.image_obj.get_bg_channel(channel)
521 bactoscoop_logger.debug("Background subtraction completed")
522
523 if cell.object_meshdata[channel]["object_contour"]:
524
525 # Extract required cellular/object properties to extract the features
526 step_length = u.get_step_length(cell.x1, cell.y1, cell.x2, cell.y2, px)
527 cell_volume = u.get_volume(
528 cell.x1, cell.y1, cell.x2, cell.y2, step_length, px
529 )
530 cell_area = u.get_area(cell.contour, px)
531 cell_widthno = u.get_width(cell.x1, cell.y1, cell.x2, cell.y2)
532 # Object properties and cropping to speed up analysis
533 object_contours = cell.object_meshdata[channel]["object_contour"]
534 object_meshes = cell.object_meshdata[channel]["object_mesh"]
535 object_midlines = cell.object_meshdata[channel]["object_midline"]
536 (
537 object_masks,
538 cropped_signal,
539 cropped_oc,
540 cropped_oml,
541 cropped_oms,
542 x_offset,
543 y_offset,
544 ) = u.get_object_masks(
545 cell.contour,
546 object_contours,
547 object_midlines,
548 object_meshes,
549 self.image_obj.bg_channels[channel],
550 )
551
552 cropped_contour, cropped_ml, cropped_msh = u.get_cropped_cell_data(
553 cell.contour, cell.midline, cell.mesh, x_offset, y_offset
554 )
555 # Calculations without object mesh
556 (
557 rectl,
558 rectw,
559 cell.objects_features["cell_total_obj_rect_length"],
560 cell.objects_features["cell_avg_obj_rect_width"],
561 ) = u.get_object_rect_length(object_contours, px)
562 aspect_r, cell.objects_features["cell_avg_obj_aspectratio"] = (
563 u.get_object_aspect_ratio(rectl, rectw)
564 )
565 cell.objects_features["cell_avg_obj_width_dt"], width_no_dt = (
566 u.get_object_avg_width(object_masks, px)
567 )
568 cell.objects_features["object_number"] = len(object_contours)
569
570 # Object area and NC ratio
571 area_list, total_obj_area = u.get_object_area(object_contours, px)
572 cell.objects_features["cell_total_obj_area"] = total_obj_area
573 cell.objects_features["NC_ratio"] = total_obj_area / cell_area
574
575 # Curvature characteristics
576 (
577 curvatures,
578 max_c,
579 min_c,
580 mean_c,
581 std_c,
582 cell.objects_features["cell_avg_obj_mean_c"],
583 cell.objects_features["cell_avg_obj_std_c"],
584 ) = u.get_object_curvature_characteristics(object_contours, px)
585
586 # Perimeter characteristics
587 (
588 perimeters,
589 circularities,
590 compactnesses,
591 sphericities,
592 avg_peri,
593 avg_cir,
594 avg_comp,
595 avg_sph,
596 ) = u.get_object_perimeter_measurements(object_contours, area_list, px)
597 cell.objects_features["cell_avg_obj_perimeter"] = avg_peri
598 cell.objects_features["cell_avg_obj_circularity"] = avg_cir
599 cell.objects_features["cell_avg_obj_compactness"] = avg_comp
600 cell.objects_features["cell_avg_obj_sphericity"] = avg_sph
601
602 (
603 convexity,
604 eccentricity,
605 solidity,
606 avg_convexity,
607 avg_eccentricity,
608 avg_solidity,
609 ) = u.get_object_convexity(object_masks)
610 cell.objects_features["cell_avg_obj_convexity"] = avg_convexity
611 cell.objects_features["cell_avg_obj_eccentricity"] = avg_eccentricity
612 cell.objects_features["cell_avg_obj_solidity"] = avg_solidity
613
614 # Object volume using approximation method (https://doi.org/10.1016/j.cell.2021.05.037)
615 approx_object_volumes, total_volume = u.get_approx_nucleoid_volume(
616 cell_volume, area_list, cell_area
617 )
618 cell.objects_features["cell_total_obj_volume_approx"] = total_volume
619
620 (
621 contour_bending_energies,
622 cell.objects_features["cell_obj_contour_avg_bending_energy"],
623 ) = u.object_bending_energy(object_contours, px)
624
625 # A bad object mesh does not invalidate its detected contour. Check
626 # every mesh before calculating any mesh-dependent aggregate.
627 mesh_data_complete = (
628 len(object_meshes) == len(object_contours)
629 and len(object_midlines) == len(object_contours)
630 and all(np.isfinite(area) and area > 0 for area in area_list)
631 and all(
632 np.asarray(mesh).ndim == 2
633 and np.asarray(mesh).shape[0] >= 3
634 and np.asarray(mesh).shape[1] == 4
635 and np.all(np.isfinite(mesh))
636 and np.any(np.asarray(mesh)[:, 0])
637 and np.asarray(midline).ndim == 2
638 and np.asarray(midline).shape[0] >= 2
639 and np.asarray(midline).shape[1] == 2
640 and np.all(np.isfinite(midline))
641 and np.linalg.norm(midline[0] - midline[-1]) > 1e-9
642 for mesh, midline in zip(object_meshes, object_midlines)
643 )
644 )
645 meshes_valid = False
646 ratio_list = []
647 if mesh_data_complete:
648 obj_x1, obj_y1, obj_x2, obj_y2 = u.get_obj_mesh_coords(cropped_oms)
649 obj_step_lengths, obj_lengths = u.get_object_step_length(
650 obj_x1, obj_y1, obj_x2, obj_y2, px
651 )
652 width_list, width_no, avg_obj_width = u.get_obj_avg_width(
653 obj_x1, obj_y1, obj_x2, obj_y2, px
654 )
655 ratio_list = u.ratio_mesh_over_total_area(
656 obj_step_lengths, width_no, area_list
657 )
658 meshes_valid = len(ratio_list) == len(object_contours) and all(
659 np.isfinite(ratio) and 0.9 <= ratio <= 1.1
660 for ratio in ratio_list
661 )
662
663 mesh_feature_names = (
664 "cell_total_obj_mesh_length",
665 "cell_avg_obj_width",
666 "cell_avg_obj_width_variability",
667 "object_avg_consDegree",
668 "cell_total_obj_volume",
669 "cell_total_obj_surface_area",
670 "cell_avg_obj_SOV",
671 "cell_obj_midline_avg_bending_energy",
672 "cell_avg_midline_sinuosity",
673 )
674 if meshes_valid:
675 cropped_signal_interp2d = u.interp2d(cropped_signal)
676 axial_intensity = u.measure_along_objects_midline_interp2d(
677 cropped_oml, cropped_signal_interp2d, width=5
678 )
679 cell.objects_features["cell_total_obj_mesh_length"] = sum(obj_lengths)
680 cell.objects_features["cell_avg_obj_width"] = avg_obj_width
681 width_var_list, cell.objects_features["cell_avg_obj_width_variability"] = (
682 u.get_object_width_variability(width_no)
683 )
684 (
685 constrDegree,
686 relPos,
687 absconstrDegree,
688 ctpos,
689 cell.objects_features["object_avg_consDegree"],
690 avg_abs_consDegree,
691 ) = u.constr_degree(axial_intensity, obj_step_lengths, width_no)
692 object_volumes, cell.objects_features["cell_total_obj_volume"] = (
693 u.get_nucleoid_volume(
694 obj_x1, obj_y1, obj_x2, obj_y2, obj_step_lengths, px
695 )
696 )
697 surface_area_list, cell.objects_features["cell_total_obj_surface_area"] = (
698 u.get_object_surface_area(width_no, obj_step_lengths)
699 )
700 sov_list, cell.objects_features["cell_avg_obj_SOV"] = (
701 u.get_object_surface_area_over_volume(
702 surface_area_list, object_volumes
703 )
704 )
705 (
706 midline_bending_energies,
707 cell.objects_features["cell_obj_midline_avg_bending_energy"],
708 ) = u.object_bending_energy(object_midlines, px)
709 sinuosity_list, cell.objects_features["cell_avg_midline_sinuosity"] = (
710 u.get_object_sinuosity(object_midlines)
711 )
712 else:
713 bactoscoop_logger.debug(
714 "Object mesh quality failed for cell %s, channel %s "
715 "(complete=%s, ratios=%s); "
716 "%s",
717 cell.cell_id, channel, mesh_data_complete, ratio_list,
718 "keeping contour features and omitting mesh features"
719 if retain_contour_on_object_mesh_failure
720 else "omitting the objects feature row",
721 )
722 if not retain_contour_on_object_mesh_failure:
723 cell.objects_features = {}
724 cell.discard = True
725 return
726 cell.objects_features.update(
727 {name: np.nan for name in mesh_feature_names}
728 )
729 sov_list = constrDegree = relPos = midline_bending_energies = sinuosity_list = np.nan
730
731 object_centroid_list = u.get_object_centroid(object_masks)
732 (
733 cell.objects_features["cell_avg_interobject_distance"],
734 interobject_distances,
735 ) = u.get_interobject_distance(object_centroid_list, px)
736
737 pole1, pole2 = u.get_pole_coordinates(cropped_msh)
738
739 (
740 cell.objects_features["avg_pole_obj_distance"],
741 min2distances,
742 pole_distances,
743 ) = u.get_pole_object_distances(cropped_msh, object_centroid_list, px)
744 (
745 cell.objects_features["l"],
746 cell.objects_features["d"],
747 cell.objects_features["l_norm"],
748 cell.objects_features["d_norm"],
749 ) = u.projection_on_midline(
750 object_centroid_list, cropped_ml, cell_widthno, px
751 )
752 cell.objects_features["object_avg_center_distance"] = (
753 u.get_avg_distance_from_center(cell.objects_features["l_norm"])
754 )
755
756 # Texture features
757 (
758 glcm_object_dissimilarity,
759 glcm_object_correlation,
760 glcm_object_homogeneity,
761 glcm_object_energy,
762 glcm_object_contrast,
763 glcm_object_entropy,
764 ) = u.calculate_object_haralick_features(object_masks, cropped_signal)
765 cell.objects_features["glcm_obj_dissimilarity"] = glcm_object_dissimilarity
766 cell.objects_features["glcm_obj_correlation"] = glcm_object_correlation
767 cell.objects_features["glcm_obj_homogeneity"] = glcm_object_homogeneity
768 cell.objects_features["glcm_obj_energy"] = glcm_object_energy
769 cell.objects_features["glcm_obj_contrast"] = glcm_object_contrast
770 cell.objects_features["glcm_obj_shannon_entropy"] = glcm_object_entropy
771
772 object_intensity_values, object_intensities_list = u.get_object_intensities(
773 object_masks, cropped_signal
774 )
775 cell.objects_features["objects_intensity_mean"] = np.mean(
776 object_intensity_values
777 )
778 cell.objects_features["objects_intensity_std"] = np.std(
779 object_intensity_values
780 )
781 cell.objects_features["objects_kurtosis"] = u.get_kurtosis(
782 object_intensity_values
783 )
784 cell.objects_features["objects_skewness"] = u.get_skew(
785 object_intensity_values
786 )
787
788 if all_data:
789 cell.objects_features["obj_areas"] = area_list
790 cell.objects_features["obj_aspect_ratios"] = aspect_r
791 cell.objects_features["max_obj_curvatures"] = max_c
792 cell.objects_features["min_obj_curvatures"] = min_c
793 cell.objects_features["mean_obj_curvatures"] = mean_c
794 cell.objects_features["std_obj_curvatures"] = std_c
795 cell.objects_features["obj_perimeters"] = perimeters
796 cell.objects_features["obj_circularities"] = circularities
797 cell.objects_features["obj_compactnesses"] = compactnesses
798 cell.objects_features["obj_sphericities"] = sphericities
799 cell.objects_features["obj_convexities"] = convexity
800 cell.objects_features["obj_eccentricities"] = eccentricity
801 cell.objects_features["obj_solidities"] = solidity
802 cell.objects_features["obj_volumes_approx"] = approx_object_volumes
803 cell.objects_features["obj_SOV"] = sov_list
804 cell.objects_features["obj_constriction_degrees"] = constrDegree
805 cell.objects_features["obj_relative_positions"] = relPos
806 cell.objects_features["obj_contour_bending_energies"] = (
807 contour_bending_energies
808 )
809 cell.objects_features["obj_midline_bending_energies"] = (
810 midline_bending_energies
811 )
812 cell.objects_features["obj_interobject_distances"] = (
813 interobject_distances
814 )
815 cell.objects_features["obj_midline_sinuosities"] = sinuosity_list
816 cell.objects_features["min_2_distances_from_pole"] = min2distances
817 cell.objects_features["obj_pole_distances"] = pole_distances
818 cell.objects_features["obj_intensity_values"] = object_intensity_values
819 cell.objects_features["obj_mean_intensities"] = u.get_objects_mean(
820 object_intensities_list
821 )
822 cell.objects_features["obj_std_intensities"] = u.get_objects_std(
823 object_intensities_list
824 )
825 cell.objects_features["obj_kurtosis_intensities"] = (
826 u.get_objects_kurtosis(object_intensities_list)
827 )
828 cell.objects_features["obj_skewness_intensities"] = (
829 u.get_objects_skewness(object_intensities_list)
830 )
831
832 else:
833 raise ValueError("No object contours to calculate features from")
834 cell.discard = True
835
836 def membrane(self, cell, channel, all_data, use_shifted_contours, shift_signal):
837 cell.membrane_features = {}
838 px = self.image_obj.px
839
840 if self.image_obj.channels is None or channel not in self.image_obj.channels:
841 raise AttributeError(
842 f"No objects are detected for channel {channel} or the channel is not loaded"
843 )
844 return
845
846 if channel not in self.image_obj.bg_channels:
847 self.image_obj.get_bg_channel(channel)
848 bactoscoop_logger.debug("Background subtraction completed")
849
850 if cell.object_meshdata[channel]["object_contour"]:
851
852 contour, midline, mesh = self.contour_selection(cell, use_shifted_contours)
853
854 if shift_signal:
855 cropped_signal, cropped_mask, _, x_offset, y_offset = u.crop_image(
856 self.image_obj.bg_channels[channel],
857 contour,
858 phase_img=self.image_obj.image,
859 )
860 else:
861 cropped_signal, cropped_mask, _, x_offset, y_offset = u.crop_image(
862 self.image_obj.bg_channels[channel],
863 contour,
864 phase_img=None
865 )
866
867 cropped_contour, _, cropped_msh = u.get_cropped_cell_data(
868 contour, midline, mesh, x_offset, y_offset
869 )
870 cropped_signal_interp2d = u.interp2d(cropped_signal)
871
872 p2p_contour = u.mesh2contour(
873 cropped_msh[:, 0], cropped_msh[:, 1], cropped_msh[:, 2], cropped_msh[:, 3]
874 )
875 eroded_contour = u.erode_contour(p2p_contour, scale=2)
876 eroded_contour_aligned = u.align_contour_start(eroded_contour, [cropped_msh[0, 0], cropped_msh[0, 1]])
877 area = u.get_area(cropped_contour, px)
878 cell.membrane_features["cell_perimeter"], _, _, _ = (
879 u.get_cell_perimeter_measurements(cropped_contour, area, px)
880 )
881 binned_eroded_contour = u.change_contour_array_length(
882 cell.membrane_features["cell_perimeter"], eroded_contour_aligned
883 )
884
885 contour_intensities = u.measure_intensity_interp2d(
886 binned_eroded_contour, cropped_signal_interp2d, width=6
887 )
888 contour_intensity_profile = u.gaussian_smoothing(contour_intensities)
889 cell.membrane_features["contour_intensity"] = contour_intensity_profile
890 cell.membrane_features["normalized_contour_intensity"] = u.normalize_per_cell(
891 cell.membrane_features["contour_intensity"]
892 )
893 cell.membrane_features["complemented_contour_intensity"] = (
894 u.get_complemented_contour_intensity(
895 cell.membrane_features["normalized_contour_intensity"]
896 )
897 )
898
899 cell.membrane_features["contour_intensity_mean"] = np.mean(
900 contour_intensity_profile
901 )
902 cell.membrane_features["contour_intensity_std"] = np.std(
903 contour_intensity_profile
904 )
905 cell.membrane_features["contour_intensity_kurtosis"] = u.get_kurtosis(
906 contour_intensity_profile
907 )
908 cell.membrane_features["contour_intensity_skewness"] = u.get_skew(
909 contour_intensity_profile
910 )
911 cell.membrane_features["contour_gradient_homogeneity"] = u.get_homogeneity(
912 cell.membrane_features["normalized_contour_intensity"]
913 )
914
915 (
916 glcm_contour_dissimilarity,
917 glcm_contour_correlation,
918 glcm_contour_homogeneity,
919 glcm_contour_energy,
920 glcm_contour_contrast,
921 glcm_contour_entropy,
922 ) = u.calculate_contour_haralick_features(contour_intensities)
923 cell.membrane_features["glcm_contour_dissimilarity"] = (
924 glcm_contour_dissimilarity
925 )
926 cell.membrane_features["glcm_contour_correlation"] = glcm_contour_correlation
927 cell.membrane_features["glcm_contour_homogeneity"] = glcm_contour_homogeneity
928 cell.membrane_features["glcm_contour_energy"] = glcm_contour_energy
929 cell.membrane_features["glcm_contour_contrast"] = glcm_contour_contrast
930 cell.membrane_features["glcm_contour_shannon_entropy"] = glcm_contour_entropy
931
932 self.discard_bad_meshes(cell)
933 else:
934 raise ValueError("No object contours to calculate features from")
935 cell.discard = True
936#not yet functional
937 def colocalization(
938 self,
939 cell,
940 signal1,
941 signal2,
942 shift_signal1,
943 shift_signal2,
944 all_data,
945 use_shifted_contours,
946 ):
947
948 cell.signal_correlation_features = {}
949 px = self.image_obj.px
950
951 if (
952 self.image_obj.channels is None
953 or signal1 not in self.image_obj.channels
954 or signal2 not in self.image_obj.channels
955 ):
956 raise AttributeError(
957 f"Channel {signal1} or {signal2} is not detected or not loaded"
958 )
959
960 if signal1 not in self.image_obj.bg_channels:
961 self.image_obj.get_bg_channel(signal1)
962 bactoscoop_logger.debug("Background subtraction completed for signal1")
963
964 if signal2 not in self.image_obj.bg_channels:
965 self.image_obj.get_bg_channel(signal2)
966 bactoscoop_logger.debug("Background subtraction completed for signal2")
967
968 contour, midline, mesh = self.contour_selection(cell, use_shifted_contours)
969
970 cropped_signal1, cropped_mask1, _, x_offset1, y_offset1 = (
971 u.crop_signal_with_optional_shift(
972 self.image_obj, signal1, contour, shift_signal1
973 )
974 )
975 cropped_signal2, cropped_mask2, _, x_offset2, y_offset2 = (
976 u.crop_signal_with_optional_shift(
977 self.image_obj, signal2, contour, shift_signal2
978 )
979 )