  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        )