|
1 | | -"""Texture featurization module scaffold.""" |
| 1 | +"""This module generates texture features for each object in the |
| 2 | +image using Haralick features. |
2 | 3 |
|
3 | | -from __future__ import annotations |
| 4 | +We do this in a as close to zero-copy way as possible. |
| 5 | +We want to make this module fast, memory efficient, and robust to large images |
| 6 | +and objects. |
| 7 | +We want this module to be python api callable and scalable. |
| 8 | +""" |
4 | 9 |
|
5 | | -from zedprofiler.exceptions import ZedProfilerError |
| 10 | +import mahotas |
| 11 | +import numpy |
| 12 | +import skimage |
| 13 | +import skimage.measure |
6 | 14 |
|
| 15 | +from zedprofiler.IO.loading_classes import ObjectLoader |
7 | 16 |
|
8 | | -def compute() -> dict[str, list[float]]: |
9 | | - """Placeholder for texture computation implementation.""" |
10 | | - raise ZedProfilerError("texture.compute is not implemented yet") |
| 17 | + |
| 18 | +def scale_image(image: numpy.ndarray, num_gray_levels: int = 256) -> numpy.ndarray: |
| 19 | + """ |
| 20 | + Scale the image to a specified number of gray levels. |
| 21 | + Example: 1024 gray levels will be scaled to 256 gray levels if |
| 22 | + num_gray_levels=256. |
| 23 | + An image with a pixel value of 0 will be scaled to 0 and a pixel value |
| 24 | + of 1023 will be scaled to 255. |
| 25 | +
|
| 26 | + Parameters |
| 27 | + ---------- |
| 28 | + image : numpy.ndarray |
| 29 | + The input image to be scaled. Can be a ndarray of any shape. |
| 30 | + num_gray_levels : int, optional |
| 31 | + The number of gray levels to scale the image to, by default 256 |
| 32 | +
|
| 33 | + Returns |
| 34 | + ------- |
| 35 | + numpy.ndarray |
| 36 | + The gray level scaled image of any shape. |
| 37 | + """ |
| 38 | + outrange_mapping = { |
| 39 | + 256: "uint8", |
| 40 | + 65536: "uint16", |
| 41 | + } |
| 42 | + try: |
| 43 | + out_range = outrange_mapping.get(num_gray_levels) |
| 44 | + except KeyError: |
| 45 | + out_range = None |
| 46 | + if out_range is None: |
| 47 | + raise ValueError( |
| 48 | + f"Unsupported num_gray_levels: {num_gray_levels}. " |
| 49 | + f"Supported values are: {list(outrange_mapping.keys())}" |
| 50 | + ) |
| 51 | + # scale the image to the requested gray levels |
| 52 | + return skimage.exposure.rescale_intensity( |
| 53 | + image, |
| 54 | + in_range="image", |
| 55 | + out_range=out_range, |
| 56 | + ) |
| 57 | + |
| 58 | + |
| 59 | +def compute_texture( |
| 60 | + object_loader: ObjectLoader, |
| 61 | + distance: int = 1, |
| 62 | + grayscale: int = 256, |
| 63 | +) -> dict: |
| 64 | + """ |
| 65 | + Calculate texture features for each object in the image using Haralick features. |
| 66 | +
|
| 67 | + The features are calculated for each object separately and the mean value |
| 68 | + is returned. |
| 69 | +
|
| 70 | + Parameters |
| 71 | + ---------- |
| 72 | + object_loader : ObjectLoader |
| 73 | + The object loader containing the image and object information. |
| 74 | + distance : int, optional |
| 75 | + The distance parameter for Haralick features, by default 1 |
| 76 | + grayscale : int, optional |
| 77 | + The number of gray levels to scale the image to, by default 256 |
| 78 | +
|
| 79 | + Returns |
| 80 | + ------- |
| 81 | + dict |
| 82 | + A dictionary containing the object ID, texture name, and texture value |
| 83 | + with keys: |
| 84 | + - object_id |
| 85 | + - texture_name |
| 86 | + - texture_value |
| 87 | +
|
| 88 | + Texture names include: Angular Second Moment, Contrast, Correlation, |
| 89 | + Variance, Inverse Difference Moment, Sum Average, Sum Variance, |
| 90 | + Sum Entropy, Entropy, and related texture measures. |
| 91 | +
|
| 92 | + - AngularSecondMoment |
| 93 | + - Contrast |
| 94 | + - Correlation |
| 95 | + - Variance |
| 96 | + - InverseDifferenceMoment |
| 97 | + - SumAverage |
| 98 | + - SumVariance |
| 99 | + - SumEntropy |
| 100 | + - Entropy |
| 101 | + - DifferenceVariance |
| 102 | + - DifferenceEntropy |
| 103 | + - InformationMeasureOfCorrelation1 |
| 104 | + - InformationMeasureOfCorrelation2 |
| 105 | +
|
| 106 | + """ |
| 107 | + label_object = object_loader.label_image |
| 108 | + labels = object_loader.object_ids |
| 109 | + feature_names = [ |
| 110 | + "AngularSecondMoment", |
| 111 | + "Contrast", |
| 112 | + "Correlation", |
| 113 | + "Variance", |
| 114 | + "InverseDifferenceMoment", |
| 115 | + "SumAverage", |
| 116 | + "SumVariance", |
| 117 | + "SumEntropy", |
| 118 | + "Entropy", |
| 119 | + "DifferenceVariance", |
| 120 | + "DifferenceEntropy", |
| 121 | + "InformationMeasureOfCorrelation1", |
| 122 | + "InformationMeasureOfCorrelation2", |
| 123 | + ] |
| 124 | + # set the number of directions based on the dimensionality of the image |
| 125 | + n_directions = 13 |
| 126 | + |
| 127 | + output_texture_dict = { |
| 128 | + "object_id": [], |
| 129 | + "texture_name": [], |
| 130 | + "texture_value": [], |
| 131 | + } |
| 132 | + # Precompute bboxes for labeled regions to avoid per-object full-array copies. |
| 133 | + props = skimage.measure.regionprops_table( |
| 134 | + label_object, |
| 135 | + properties=["label", "bbox"], |
| 136 | + ) |
| 137 | + # Map label id to bbox (z0, y0, x0, z1, y1, x1) |
| 138 | + label_to_bbox = {} |
| 139 | + labels_prop = props.get("label", []) |
| 140 | + for i, lbl in enumerate(labels_prop): |
| 141 | + label_to_bbox[int(lbl)] = ( |
| 142 | + int(props["bbox-0"][i]), |
| 143 | + int(props["bbox-1"][i]), |
| 144 | + int(props["bbox-2"][i]), |
| 145 | + int(props["bbox-3"][i]), |
| 146 | + int(props["bbox-4"][i]), |
| 147 | + int(props["bbox-5"][i]), |
| 148 | + ) |
| 149 | + # loop through each label and get the bounding box |
| 150 | + # to compute features for the object |
| 151 | + for _, label in enumerate(labels): |
| 152 | + if int(label) == 0: |
| 153 | + continue |
| 154 | + bbox = label_to_bbox.get(int(label)) |
| 155 | + if bbox is None: |
| 156 | + continue |
| 157 | + |
| 158 | + min_z, min_y, min_x, max_z, max_y, max_x = bbox |
| 159 | + |
| 160 | + # Crop to the object's bounding box (skimage bboxes are half-open) |
| 161 | + image_object = object_loader.image[min_z:max_z, min_y:max_y, min_x:max_x].copy() |
| 162 | + selected_label_object = label_object[min_z:max_z, min_y:max_y, min_x:max_x] |
| 163 | + object_mask = selected_label_object == label |
| 164 | + if not numpy.any(object_mask): |
| 165 | + continue |
| 166 | + image_object[~object_mask] = 0 |
| 167 | + features = numpy.empty((n_directions, 13, max(labels))) |
| 168 | + image_object = scale_image(image_object, num_gray_levels=grayscale) |
| 169 | + try: |
| 170 | + # calculates 13 Haralick features for each direction (13) |
| 171 | + # and each object, and stores them in a 3D array |
| 172 | + features[:, :, label - 1] = mahotas.features.haralick( |
| 173 | + ignore_zeros=True, |
| 174 | + f=image_object, |
| 175 | + distance=distance, |
| 176 | + compute_14th_feature=False, |
| 177 | + ) |
| 178 | + except ValueError: |
| 179 | + features = numpy.full(len(feature_names), numpy.nan, dtype=float) |
| 180 | + # iterate through the direction, feature, and object dimensions |
| 181 | + # of the features array to populate the output dictionary |
| 182 | + for direction, direction_features in enumerate(features): |
| 183 | + direction_str = f"{direction:02d}" |
| 184 | + for feature_name, feature in zip(feature_names, direction_features): |
| 185 | + for object_id, feature_value in zip(labels, feature): |
| 186 | + output_texture_dict["object_id"].append(object_id) |
| 187 | + output_texture_dict["texture_name"].append( |
| 188 | + f"{feature_name}-{distance}-{direction_str}-{grayscale}" |
| 189 | + ) |
| 190 | + output_texture_dict["texture_value"].append(feature_value) |
| 191 | + return output_texture_dict |
0 commit comments