Coverage for src/recon3d/instance_analysis.py: 85%
239 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-02 00:06 +0000
1"""
2Basic workflow for processing grayscale image stacks to labeled images for
3feature property measurement.
5This module supports multiple phases (>2 phases) and includes simple
6segmentation for demonstration purposes.
8It writes outputs as image directories, but can be integrated with HDF for
9internal dataset storage.
10"""
12# TODO Integrate with HDF writing
13# TODO Integrate with feature measurements
14# TODO Use constants for Z axis in other module
15# TODO Simplify/Impose structure for image directory (dataset directory in HDF)
16# to remove magic numbers and indexing of phase_paths and label_paths
17# TODO Devise better solution for multiphase support without if/else blocks
20import argparse
21import glob
22from typing import Final, NamedTuple, Tuple, Type, List
23import math
25import h5py
26import numpy as np
27import numpy.typing as npt
28from pathlib import Path
29from PIL import Image
30import skimage.io as skio
31from sklearn.neighbors import NearestNeighbors as sklearnNN
33from scipy import ndimage
34from skimage.measure import label, regionprops
35import cc3d
36import fastremap
38import recon3d.constants as cs
39import recon3d.types
40import recon3d.utility as ut
41import recon3d.hdf_io as hio
42from recon3d.types import *
45# # ----------------------------------
46# # USER INPUT DATA HARD CODES - begin
47# # ----------------------------------
48# # NEED TO PUT INTO A .yml file instead of hard codes
49# # input_dir: str = "/Users/apolon/Desktop/test_recon/input/"
50# # save_dir: str = "/Users/apolon/Desktop/test_recon/output/"
52# input_dir: Final[str] = "/Users/chovey/temp/test_recon/input/"
53# save_dir: Final[str] = "/Users/chovey/temp/test_recon/output/"
55# file_type: Final[str] = ".tif"
56# gray_path: Final[str] = f"{save_dir}/gray/"
57# threshold: Final[int] = 185
58# segmented_path: Final[str] = f"{save_dir}/segmented/"
59# # first phase is always assumed to be the sample (1 or 255 in segmented images)
60# phase_paths: np.ndarray = np.array(
61# [f"{save_dir}/sample_mask/", f"{save_dir}/void_mask/"]
62# )
63# # first label path is always assumed to the first phase after the sample,
64# # same ordering as phase_paths
65# label_paths: np.ndarray = np.array([f"{save_dir}/void_ids/"])
67# # --------------------------------
68# # USER INPUT DATA HARD CODES - end
69# # --------------------------------
72def calc_moment(
73 indices: np.ndarray[np.int32],
74 p: int,
75 q: int,
76 r: int,
77 centroid_px: Centroid,
78) -> float:
79 """
80 Calculate 3D moment invariants of a feature based on its voxel locations.
81 Mask is given in index locations, so math here is all done in pixel units.
83 Parameters
84 ----------
85 indices : np.ndarray[np.int32]
86 Array of [Z,Y,X] indices corresponding to the feature's voxel
87 locations.
88 The shape of the array should be (n, 3), where n is the number of
89 voxels of the feature.
90 Example: np.array([[Z1, Y1, X1], [Z2, Y2, X2], ..., [Zn, Yn, Xn]]).
91 p : int
92 Order of the moment for the X dimension.
93 q : int
94 Order of the moment for the Y dimension.
95 r : int
96 Order of the moment for the Z dimension.
97 centroid_px : Centroid
98 Centroid location for X, Y, and Z axis, respectively, in pixel units.
100 Returns
101 -------
102 float
103 The 3D central moment (translation invariant) in pixel units
104 normalized by object volume.
106 Raises
107 ------
108 ValueError
109 If the units of the centroid are not in pixel/voxel units.
111 Examples
112 --------
113 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32)
114 >>> centroid_px = Centroid(
115 ... cx=Length(1.0, Units.VOXEL),
116 ... cy=Length(1.0, Units.VOXEL),
117 ... cz=Length(1.0, Units.VOXEL)
118 ... )
119 >>> calc_moment(indices, 2, 2, 2, centroid_px)
120 0.66
121 """
123 unit = centroid_px.cx.unit
124 if unit != Units.VOXEL:
125 raise ValueError(
126 f"Moment invariants are calculated using pixel/voxel units, but {unit} units were given"
127 )
129 mu = 0
130 for i in range(0, indices.shape[0]):
131 mu += (
132 (indices[i][2] - centroid_px.cx.value)
133 ** p # x-location of voxel w.r.t x-centroid
134 * (indices[i][1] - centroid_px.cy.value)
135 ** q # y-location of voxel w.r.t y-centroid
136 * (indices[i][0] - centroid_px.cz.value)
137 ** r # z-location of voxel w.r.t z-centroid
138 )
139 return mu / len(indices)
142def center_of_mass(
143 indices: np.ndarray[np.int32], resolution: Resolution, origin: Origin
144) -> Centroid:
145 """
146 Calculates the center of mass (centroid) of an instance object.
147 Requires that the instance object has uniform density (is homogeneous).
149 Parameters
150 ----------
151 indices : np.ndarray[np.int32]
152 Array of [Z,Y,X] indices corresponding to the feature's voxel
153 locations.
154 The shape of the array should be (n, 3), where n is the number of
155 voxels of the feature.
156 Example: np.array([[Z1, Y1, X1], [Z2, Y2, X2], ..., [Zn, Yn, Xn]]).
157 resolution : Resolution
158 Resolution of the voxel grid in each dimension (dx, dy, dz).
159 origin : Origin
160 Origin of the voxel grid in each dimension (x0, y0, z0).
162 Returns
163 -------
164 Centroid
165 The centroid location for X, Y, and Z axis, respectively, in microns.
167 Examples
168 --------
169 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32)
170 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON))
171 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON))
172 >>> center_of_mass(indices, resolution, origin)
173 Centroid(cx=Length(value=1.0, unit=Units.MICRON), cy=Length(value=1.0, unit=Units.MICRON), cz=Length(value=1.0, unit=Units.MICRON))
174 """
176 centroid_zyx = np.add(
177 np.multiply(
178 np.mean(indices, axis=0),
179 np.array(
180 [
181 resolution.dz.value,
182 resolution.dy.value,
183 resolution.dx.value,
184 ]
185 ),
186 ),
187 np.array(
188 [
189 origin.z0.value,
190 origin.y0.value,
191 origin.x0.value,
192 ]
193 ),
194 )
195 cx = centroid_zyx[2]
196 cy = centroid_zyx[1]
197 cz = centroid_zyx[0]
198 unit = resolution.dx.unit
199 return Centroid(
200 cx=Length(cx, unit=unit),
201 cy=Length(cy, unit=unit),
202 cz=Length(cz, unit=unit),
203 )
206def centroid_pix(
207 centroid: Centroid, resolution: Resolution, origin: Origin
208) -> Centroid:
209 """
210 Convert between real space units and pixel units.
212 Parameters
213 ----------
214 centroid : Centroid
215 The centroid location in real space units.
216 resolution : Resolution
217 Resolution of the voxel grid in each dimension (dx, dy, dz).
218 origin : Origin
219 Origin of the voxel grid in each dimension (x0, y0, z0).
221 Returns
222 -------
223 Centroid
224 The centroid location in pixel units.
226 Examples
227 --------
228 >>> centroid = Centroid(cx=Length(10.0, Units.MICRON), cy=Length(20.0, Units.MICRON), cz=Length(30.0, Units.MICRON))
229 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON))
230 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON))
231 >>> centroid_pix(centroid, resolution, origin)
232 Centroid(cx=Length(value=10.0, unit=Units.VOXEL), cy=Length(value=20.0, unit=Units.VOXEL), cz=Length(value=30.0, unit=Units.VOXEL))
233 """
235 # convert centroids to normalized pixel units
236 cx = (centroid.cx.value - origin.x0.value) / resolution.dx.value
237 cy = (centroid.cy.value - origin.y0.value) / resolution.dy.value
238 cz = (centroid.cz.value - origin.z0.value) / resolution.dz.value
240 unit = Units.VOXEL
242 return Centroid(
243 cx=Length(cx, unit=unit),
244 cy=Length(cy, unit=unit),
245 cz=Length(cz, unit=unit),
246 )
249def ellipsoid_surface_area(ellipsoid: BestFitEllipsoid) -> Area:
250 """
251 Calculate the surface area of the best fit ellipsoid using the Knud Thomsen formula.
252 Applicable to scalene ellipsoids (semi-axes a > b > c), with a relative error of at most 1.061%.
254 Parameters
255 ----------
256 ellipsoid : BestFitEllipsoid
257 The best fit ellipsoid with semi-axes lengths a, b, and c.
259 Returns
260 -------
261 Area
262 The surface area of the ellipsoid in microns squared.
264 Examples
265 --------
266 >>> ellipsoid = BestFitEllipsoid(
267 ... a=EllipsoidAxis(length=Length(5.0, Units.MICRON), orientation=UnitVector(u=1, v=0, w=0)),
268 ... b=EllipsoidAxis(length=Length(3.0, Units.MICRON), orientation=UnitVector(u=0, v=1, w=0)),
269 ... c=EllipsoidAxis(length=Length(2.0, Units.MICRON), orientation=UnitVector(u=0, v=0, w=1))
270 ... )
271 >>> ellipsoid_surface_area(ellipsoid)
272 Area(value=134.8149867941625, unit_squared=<Units.MICRON: 'micron'>)
273 """
274 unit = ellipsoid.a.length.unit
275 p = cs.Constants().KNUD_THOMSEN_FACTOR
276 a, b, c = (
277 ellipsoid.a.length.value,
278 ellipsoid.b.length.value,
279 ellipsoid.c.length.value,
280 )
282 surface_area = (
283 4 * math.pi * ((((a**p * b**p) + (a**p * c**p) + (b**p * c**p)) / 3) ** (1 / p))
284 )
285 return Area(value=surface_area, unit_squared=unit)
288def ellipsoid_volume(ellipsoid: BestFitEllipsoid) -> Volume:
289 """
290 Calculate the volume of the best fit ellipsoid using semi-axes lengths.
292 Parameters
293 ----------
294 ellipsoid : BestFitEllipsoid
295 The best fit ellipsoid with semi-axes lengths a, b, and c.
297 Returns
298 -------
299 Volume
300 The volume of the ellipsoid in microns cubed.
302 Examples
303 --------
304 >>> ellipsoid = BestFitEllipsoid(
305 ... a=EllipsoidAxis(length=Length(5.0, Units.MICRON), orientation=UnitVector(u=1, v=0, w=0)),
306 ... b=EllipsoidAxis(length=Length(3.0, Units.MICRON), orientation=UnitVector(u=0, v=1, w=0)),
307 ... c=EllipsoidAxis(length=Length(2.0, Units.MICRON), orientation=UnitVector(u=0, v=0, w=1))
308 ... )
309 >>> ellipsoid_volume(ellipsoid)
310 Volume(value=125.66370614359171, unit_cubed=<Units.MICRON: 'micron'>)
311 """
313 unit = ellipsoid.a.length.unit
314 a, b, c = (
315 ellipsoid.a.length.value,
316 ellipsoid.b.length.value,
317 ellipsoid.c.length.value,
318 )
320 volume = 4 / 3 * math.pi * a * b * c
322 return Volume(value=volume, unit_cubed=unit)
325def equivalent_spherical_diameter(
326 n_voxels: int,
327 resolution: Resolution,
328) -> Length:
329 """
330 Calculate the equivalent spherical diameter of a feature.
332 Parameters
333 ----------
334 n_voxels : int
335 The number of voxels in the feature.
336 resolution : Resolution
337 The resolution of the voxel grid in each dimension (dx, dy, dz).
339 Returns
340 -------
341 Length
342 The equivalent spherical diameter of the feature in the same units as
343 the resolution.
345 Examples
346 --------
347 >>> resolution = Resolution(
348 ... dx=Length(1.0, Units.MICRON),
349 ... dy=Length(1.0, Units.MICRON),
350 ... dz=Length(1.0, Units.MICRON)
351 ... )
352 >>> equivalent_spherical_diameter(1000, resolution)
353 Length(value=12.407009817988, unit=<Units.MICRON: 'micron'>)
354 """
356 value = 2 * np.power(
357 (
358 (3 / (4 * np.pi))
359 * (
360 n_voxels
361 * resolution.dx.value
362 * resolution.dy.value
363 * resolution.dz.value
364 )
365 ),
366 (1 / 3),
367 )
368 return Length(value=value, unit=resolution.dx.unit)
371def fit_ellipsoid(
372 indices: np.ndarray,
373 centroid: Centroid,
374 resolution: Resolution,
375 origin: Origin,
376) -> BestFitEllipsoid:
377 """
378 Fit a 3D ellipsoid to a set of voxel indices.
380 Parameters
381 ----------
382 indices : np.ndarray
383 Array of [Z,Y,X] indices corresponding to the feature's voxel
384 locations.
385 The shape of the array should be (n, 3), where n is the number of
386 voxels of the feature.
387 centroid : Centroid
388 Centroid location for X, Y, and Z axis, respectively, in spatial
389 coordinates (e.g., microns).
390 resolution : Resolution
391 Voxel resolution for X, Y, and Z dimensions, respectively, in length
392 units (e.g., microns).
393 origin : Origin
394 Optional origin in spatial coordinates if the volume has been
395 previously transformed.
397 Returns
398 -------
399 BestFitEllipsoid
400 The best fit ellipsoid with semi-axes lengths and orientations.
402 Examples
403 --------
404 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32)
405 >>> centroid = Centroid(cx=Length(1.0, Units.MICRON), cy=Length(1.0, Units.MICRON), cz=Length(1.0, Units.MICRON))
406 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON))
407 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON))
408 >>> fit_ellipsoid(indices, centroid, resolution, origin)
409 BestFitEllipsoid(a=EllipsoidAxis(length=Length(value=80.146386085398, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=-0.25175905655605574, v=-0.5305735687527386, w=-0.8093880809494217)), b=EllipsoidAxis(length=Length(value=1.2477168951003976, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=-0.8774683531474325, v=-0.227651095921704, w=0.4221661613042665)), c=EllipsoidAxis(length=Length(value=2.380014006104306e-07, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.4082482904637247, v=-0.816496580927762, w=0.4082482904639296)))
410 """
412 if not resolution.dx == resolution.dy == resolution.dz:
413 raise NotImplementedError(
414 "Only isotropic voxels are currently supported for axis length calculation"
415 )
417 unit = resolution.dx.unit
419 # convert centroids to normalized pixel units
420 cen_px = centroid_pix(centroid, resolution, origin)
422 # Get moments
423 s_xx = calc_moment(indices, 2, 0, 0, cen_px)
424 s_yy = calc_moment(indices, 0, 2, 0, cen_px)
425 s_zz = calc_moment(indices, 0, 0, 2, cen_px)
426 s_xy = calc_moment(indices, 1, 1, 0, cen_px)
427 s_xz = calc_moment(indices, 1, 0, 1, cen_px)
428 s_yz = calc_moment(indices, 0, 1, 1, cen_px)
430 # Create matrix of 3D moments
431 s_ij = np.array(
432 [
433 [s_xx, s_xy, s_xz],
434 [s_xy, s_yy, s_yz],
435 [s_xz, s_yz, s_zz],
436 ]
437 )
438 eig_val, eig_vec = np.linalg.eig(s_ij)
440 semi_ax_lengths_px = np.sqrt(5 * eig_val) # still in pixel units
441 sort = np.argsort(semi_ax_lengths_px)[::-1] # sort descending axis lengths
443 # TODO: math to account for change in length as function
444 # as function of axis orientation in anisotropic voxels
445 # remove NotImplementedError above
446 semi_ax_lengths = np.array(
447 [
448 semi_ax_lengths_px[sort[0]] * resolution.dx.value,
449 semi_ax_lengths_px[sort[1]] * resolution.dx.value,
450 semi_ax_lengths_px[sort[2]] * resolution.dx.value,
451 ]
452 )
453 # some axes lengths are NAN, should overwrite to 0.0
454 semi_ax_lengths = np.nan_to_num(semi_ax_lengths, nan=0.0)
456 # from docs on np.linalg.eig: The normalized (unit “length”) eigenvectors,
457 # such that the column eigenvectors[:,i] is the eigenvector corresponding to the eigenvalue eigenvalues[i].
458 ax_vectors = np.array(
459 [
460 eig_vec[:, sort[0]],
461 eig_vec[:, sort[1]],
462 eig_vec[:, sort[2]],
463 ]
464 )
465 ax_vectors = np.reshape(ax_vectors, (9))
467 a_vector = UnitVector(
468 u=ax_vectors[0],
469 v=ax_vectors[1],
470 w=ax_vectors[2],
471 )
472 a = EllipsoidAxis(
473 length=Length(semi_ax_lengths[0], unit=unit),
474 orientation=a_vector,
475 )
476 b_vector = UnitVector(
477 u=ax_vectors[3],
478 v=ax_vectors[4],
479 w=ax_vectors[5],
480 )
481 b = EllipsoidAxis(
482 length=Length(semi_ax_lengths[1], unit=unit),
483 orientation=b_vector,
484 )
485 c_vector = UnitVector(
486 u=ax_vectors[6],
487 v=ax_vectors[7],
488 w=ax_vectors[8],
489 )
490 c = EllipsoidAxis(
491 length=Length(semi_ax_lengths[2], unit=unit),
492 orientation=c_vector,
493 )
495 ellipsoid = BestFitEllipsoid(a=a, b=b, c=c)
496 return ellipsoid
499def instance_indices(
500 instance_stack: InstanceImageStack,
501) -> InstanceIndices:
502 """
503 Returns an array of instances and a list of N-dimensional arrays of indices
504 of the unique values in the instance stack.
506 Parameters
507 ----------
508 instance_stack : InstanceImageStack
509 Array with arbitrary dimensions.
511 Returns
512 -------
513 InstanceIndices
514 - 1D-array of sorted unique values.
515 - List of arrays. Each array contains the indices where a given value in x is found.
517 Examples
518 --------
519 >>> instance_stack = InstanceImageStack(
520 ... name="example_stack",
521 ... metadata=MetaData(
522 ... data_volume=DataVolume(x_width=2, y_height=2, z_image_count=1),
523 ... resolution=Resolution(
524 ... dx=Length(1.0, Units.MICRON),
525 ... dy=Length(1.0, Units.MICRON),
526 ... dz=Length(1.0, Units.MICRON)
527 ... ),
528 ... pixel_units=Units.MICRON,
529 ... origin=Origin(
530 ... x0=Length(0.0, Units.MICRON),
531 ... y0=Length(0.0, Units.MICRON),
532 ... z0=Length(0.0, Units.MICRON)
533 ... )
534 ... ),
535 ... data=np.array([[1, 2, 2], [3, 1, 1]]),
536 ... nlabels=3,
537 ... min_feature_size=1
538 ... )
539 >>> instance_indices(instance_stack)
540 InstanceIndices(
541 source_name='example_stack',
542 indices=[array([[0, 0], [1, 1], [1, 2]]), array([[0, 1], [0, 2]]), array([[1, 0]])],
543 labels=InstanceLabels(data=np.array([1, 2, 3]))
544 )
545 """
547 x = instance_stack.data
548 x_flat = x.ravel()
549 ix_flat = np.argsort(x_flat)
550 u, ix_u = fastremap.unique(x_flat[ix_flat], return_index=True)
551 ix_ndim = np.unravel_index(ix_flat, x.shape)
552 ix_ndim = np.c_[ix_ndim] if x.ndim > 1 else ix_flat
553 indices = np.split(ix_ndim, ix_u[1:])
554 labels = u
555 instance_ids = InstanceIndices(
556 source_name=instance_stack.name,
557 indices=indices,
558 labels=InstanceLabels(data=labels),
559 )
560 return instance_ids
563def minimum_size_filter(
564 initial_stack: InstanceImageStack,
565 initial_indices: InstanceIndices,
566) -> Tuple[InstanceImageStack, InstanceIndices]:
567 """
568 Removes features below minimum size from InstanceIndices and corresponding
569 InstanceImageStack. Should be done before shape metrics are determined.
571 Parameters
572 ----------
573 initial_stack : InstanceImageStack
574 The initial instance image stack.
575 initial_indices : InstanceIndices
576 The initial instance indices.
578 Returns
579 -------
580 Tuple[InstanceImageStack, InstanceIndices]
581 The filtered instance image stack and the filtered instance indices.
583 Raises
584 ------
585 ValueError
586 If the names of the initial stack and indices do not match.
588 Examples
589 --------
590 >>> initial_stack = InstanceImageStack(
591 ... name="example_stack",
592 ... metadata=MetaData(
593 ... data_volume=DataVolume(x_width=4, y_height=2, z_image_count=1),
594 ... resolution=Resolution(
595 ... dx=Length(1.0, Units.MICRON),
596 ... dy=Length(1.0, Units.MICRON),
597 ... dz=Length(1.0, Units.MICRON)
598 ... ),
599 ... pixel_units=Units.MICRON,
600 ... origin=Origin(
601 ... x0=Length(0.0, Units.MICRON),
602 ... y0=Length(0.0, Units.MICRON),
603 ... z0=Length(0.0, Units.MICRON)
604 ... )
605 ... ),
606 ... data=np.array([[1, 2, 2, 1], [3, 1, 1, 1]]),
607 ... nlabels=3,
608 ... min_feature_size=2
609 ... )
610 >>> initial_indices = instance_indices(initial_stack)
611 >>> minimum_size_filter(initial_stack, initial_indices)
612 2 connected components remaining after minimum size filter of 2 voxels
613 (InstanceImageStack(name='example_stack', metadata=MetaData(data_volume=DataVolume(x_width=4, y_height=2, z_image_count=1), resolution=Resolution(dx=Length(value=1.0, unit=<Units.MICRON: 'micron'>), dy=Length(value=1.0, unit=<Units.MICRON: 'micron'>), dz=Length(value=1.0, unit=<Units.MICRON: 'micron'>)), pixel_units=<Units.MICRON: 'micron'>, origin=Origin(x0=Length(value=0.0, unit=<Units.MICRON: 'micron'>), y0=Length(value=0.0, unit=<Units.MICRON: 'micron'>), z0=Length(value=0.0, unit=<Units.MICRON: 'micron'>))), data=array([[1, 0, 0, 1],[2, 1, 1, 1]], dtype=int8), nlabels=2, min_feature_size=2), InstanceIndices(source_name='example_stack', labels=InstanceLabels(data=array([0, 1, 2], dtype=int8)), indices=[array([[0, 1],[0, 2]]), array([[0, 0],[0, 3],[1, 1],[1, 2],[1, 3]]), array([[1, 0]])]))
614 """
616 if not initial_stack.name == initial_indices.source_name:
617 raise ValueError(
618 f"""Initial InstanceImageStack name ({initial_stack.name}) and
619 InstanceIndices name ({initial_indices.source_name}) do not match.
620 Function inputs must be from same dataset."""
621 )
623 n_voxels = num_voxels(initial_indices)
624 n_voxels_vals = np.asarray([i.value for i in n_voxels])
625 invalid_instance_ids = np.nonzero(n_voxels_vals < initial_stack.min_feature_size)[0]
626 masked_instance_ids = fastremap.mask(initial_stack.data, invalid_instance_ids)
627 remapped_instance_ids, _ = fastremap.renumber(
628 masked_instance_ids,
629 preserve_zero=True,
630 )
631 nlabels = np.max(np.unique(remapped_instance_ids))
633 filtered_stack = InstanceImageStack(
634 name=initial_stack.name,
635 metadata=initial_stack.metadata,
636 data=remapped_instance_ids,
637 nlabels=nlabels,
638 min_feature_size=initial_stack.min_feature_size,
639 )
641 filtered_indices = instance_indices(instance_stack=filtered_stack)
643 print(
644 f"\t\t{nlabels} connected components remaining after minimum size filter of {initial_stack.min_feature_size} voxels"
645 )
646 return filtered_stack, filtered_indices
649# TODO
650def map_features_to_voxels():
651 """To come."""
654def nearest_neighbor_distance(
655 centroids: Centroids,
656) -> NthNearestNeighbors:
657 """
658 Calculate the nearest neighbor distance for a set of centroids.
659 With nth_nearest_neighbor=1 returning the features themselves
660 (distances=0, neighbor is themselves).
662 Parameters
663 ----------
664 centroids : Centroids
665 The centroids for which to calculate the nearest neighbor distance.
667 Returns
668 -------
669 NthNearestNeighbors
670 The nth nearest neighbors, distances, and instance IDs.
672 Examples
673 --------
674 >>> centroids = Centroids(data=[
675 ... Centroid(cx=Length(0.0, Units.MICRON), cy=Length(0.0, Units.MICRON), cz=Length(0.0, Units.MICRON)),
676 ... Centroid(cx=Length(1.0, Units.MICRON), cy=Length(1.0, Units.MICRON), cz=Length(1.0, Units.MICRON)),
677 ... Centroid(cx=Length(2.0, Units.MICRON), cy=Length(2.0, Units.MICRON), cz=Length(2.0, Units.MICRON))
678 ... ])
679 >>> nearest_neighbor_distance(centroids)
680 NthNearestNeighbors(nth_nearest=2, distances=[Length(value=1.7320508075688772, unit=<Units.MICRON: 'micron'>), Length(value=1.7320508075688772, unit=<Units.MICRON: 'micron'>), Length(value=1.7320508075688772, unit=<Units.MICRON: 'micron'>)], instance_id=array([0, 0, 1]))
681 """
683 neighbor = sklearnNN(n_neighbors=2)
684 centroid_array = ut.centroids_to_ndarray(centroids=centroids)
685 neighbor.fit(centroid_array)
686 neighbor_dist, neighbor_index = neighbor.kneighbors(
687 centroid_array,
688 return_distance=True,
689 )
690 distances = neighbor_dist[:, 1]
691 instance_id = neighbor_index[:, 1]
692 distances[0], instance_id[0] = 0, 0 # set feature 0 to 0 distance and neighbor 0
694 unit = centroids.data[0].cx.unit
695 distances_list = [Length(value=i, unit=unit) for i in distances.tolist()]
696 nth_nearest = NthNearestNeighbors(
697 nth_nearest=2,
698 distances=distances_list,
699 instance_id=instance_id,
700 )
701 return nth_nearest
704def num_voxels(indices: InstanceIndices) -> list[NVoxel]:
705 """
706 Returns the number of voxels for each instance ID from an InstanceImageStack.
708 Parameters
709 ----------
710 indices : InstanceIndices
711 The instance indices containing the indices of each instance.
713 Returns
714 -------
715 List[NVoxel]
716 A list of NVoxel objects, each representing the number of voxels for an instance.
718 Examples
719 --------
720 >>> indices = InstanceIndices(
721 ... source_name="example_stack",
722 ... indices=[np.array([[0, 0], [1, 1]]), np.array([[0, 1], [1, 2]])],
723 ... labels=InstanceLabels(data=np.array([1, 2]))
724 ... )
725 >>> num_voxels(indices)
726 [NVoxel(value=2, unit=<Units.VOXEL: 'voxel'>), NVoxel(value=2, unit=<Units.VOXEL: 'voxel'>)]
727 """
729 indices = indices.indices
730 data = [NVoxel(value=len(i)) for i in indices]
731 # n_voxels = np.asarray(data, dtype=NVoxel)
732 n_voxels = data
733 return n_voxels
736def process_image_stack(yml_input_file: Path) -> SemanticImageStack:
737 """
738 Given a yml input file, converts it to a SemanticImageStack.
740 Parameters
741 ----------
742 yml_input_file : Path
743 The fully pathed yml file that specifies settings.
745 Returns
746 -------
747 SemanticImageStack
748 The SemanticImageStack object.
750 Examples
751 --------
752 >>> yml_input_file = Path("path/to/input.yml")
753 >>> semantic_image_stack = process_image_stack(yml_input_file)
754 """
756 print(f"Processing specification file: {yml_input_file}")
757 db = ut.yaml_to_dict(yml_input_file)
759 file_type = db["semantic_images_type"]
760 valid_extensions = (".tif", ".tiff")
762 if file_type not in valid_extensions:
763 raise ValueError("Error: 'image_type' must be '.tif' or '.tiff'.")
765 path_input = Path(db["semantic_images_dir"]).expanduser()
766 if not path_input.is_dir():
767 raise ValueError(f"Error, 'image_dir' of {path_input} not found.")
768 print(f"Input path: {path_input}")
770 path_output = Path(db["out_dir"]).expanduser()
771 if not path_output.is_dir():
772 Path(path_output).mkdir(
773 parents=True, exist_ok=True
774 ) # makes directory if it doesn't exist
775 print(f"Output path: {path_output}")
777 # semantic_seg_dict = db["class_labels"]
779 # # load images
780 # imageStack = skio.imread_collection(tiff_files, conserve_memory=True)
781 data = ut.read_images(path_input, file_type)
782 semantic_seg_stack_name = db["semantic_imagestack_name"]
784 # create meta data
785 (nz, ny, nx) = data.shape
786 data_volume = DataVolume(z_image_count=nz, y_height=ny, x_width=nx)
787 print(
788 f"image stack, {semantic_seg_stack_name}, has dimensions (num_images, row, col): {data.shape}"
789 )
791 pixel_units = db["pixel_units"]
792 valid_units = set(item.value for item in Units)
793 if pixel_units not in valid_units:
794 raise ValueError(
795 f"Error, '{pixel_units}' is not a valid unit, accepted units are: {valid_units}"
796 )
797 pixel_units = Units(db["pixel_units"])
799 dx, dy, dz = (
800 db["voxel_size"]["dx"],
801 db["voxel_size"]["dy"],
802 db["voxel_size"]["dz"],
803 )
804 resolution = Resolution(
805 dx=Length(dx, unit=pixel_units),
806 dy=Length(dy, unit=pixel_units),
807 dz=Length(dz, unit=pixel_units),
808 )
810 x0, y0, z0 = (
811 db["origin"]["x0"],
812 db["origin"]["y0"],
813 db["origin"]["z0"],
814 )
815 origin = Origin(
816 x0=Length(x0, unit=pixel_units),
817 y0=Length(y0, unit=pixel_units),
818 z0=Length(z0, unit=pixel_units),
819 )
821 meta = MetaData(
822 data_volume=data_volume,
823 resolution=resolution,
824 pixel_units=pixel_units,
825 origin=origin,
826 )
828 # # Optional, grey image stack
829 # # TODO: test functionality
830 # if db["grey_image_dir"]:
831 # path_grey_image_dir = Path(db["path_grey_image_dir"]).expanduser()
832 # assert path_grey_image_dir.is_dir(), "Error, 'image_dir' not found."
833 # grey_image_name = "Raw greyscale"
834 # grey_data = ut.read_images(path_grey_image_dir, file_type)
835 # assert data.shape == grey_data.shape
836 # semantic_seg_image_stack = ImageStack(
837 # name=grey_image_name, metadata=meta, data=grey_data
838 # )
840 # create SemanticImageStack
841 semantic_seg_image_stack = SemanticImageStack(
842 name=semantic_seg_stack_name, metadata=meta, data=data
843 )
845 return semantic_seg_image_stack # , path_file_output)
848def save_instance_images(
849 input_volumes: list[InstanceImageStack], save_path: Path
850) -> bool:
851 """
852 Save InstanceImageStack objects to image files.
854 This function takes a list of InstanceImageStack objects and saves each one as a series of image files.
855 The images are saved in the specified directory, with each stack saved in its own subdirectory.
857 Parameters
858 ----------
859 input_volumes : list of InstanceImageStack
860 List of InstanceImageStack objects to be saved as images.
861 save_path : Path
862 The directory where the images will be saved.
864 Returns
865 -------
866 bool
867 True if the images were saved successfully, False otherwise.
869 Examples
870 --------
871 >>> input_volumes = [
872 ... InstanceImageStack(
873 ... name="example_stack_1",
874 ... metadata=MetaData(
875 ... data_volume=DataVolume(x_width=256, y_height=256, z_image_count=10),
876 ... resolution=Resolution(
877 ... dx=Length(1.0, Units.MICRON),
878 ... dy=Length(1.0, Units.MICRON),
879 ... dz=Length(1.0, Units.MICRON)
880 ... ),
881 ... pixel_units=Units.MICRON,
882 ... origin=Origin(
883 ... x0=Length(0.0, Units.MICRON),
884 ... y0=Length(0.0, Units.MICRON),
885 ... z0=Length(0.0, Units.MICRON)
886 ... )
887 ... ),
888 ... data=np.random.rand(10, 256, 256),
889 ... nlabels=3,
890 ... min_feature_size=2
891 ... ),
892 ... InstanceImageStack(
893 ... name="example_stack_2",
894 ... metadata=MetaData(
895 ... data_volume=DataVolume(x_width=256, y_height=256, z_image_count=10),
896 ... resolution=Resolution(
897 ... dx=Length(1.0, Units.MICRON),
898 ... dy=Length(1.0, Units.MICRON),
899 ... dz=Length(1.0, Units.MICRON)
900 ... ),
901 ... pixel_units=Units.MICRON,
902 ... origin=Origin(
903 ... x0=Length(0.0, Units.MICRON),
904 ... y0=Length(0.0, Units.MICRON),
905 ... z0=Length(0.0, Units.MICRON)
906 ... )
907 ... ),
908 ... data=np.random.rand(10, 256, 256),
909 ... nlabels=3,
910 ... min_feature_size=2
911 ... )
912 ... ]
913 >>> save_path = Path("path/to/save")
914 >>> save_instance_images(input_volumes, save_path)
915 True
916 """
918 for instance_volumes in input_volumes:
919 ut.ndarray_to_img(
920 data=instance_volumes.data,
921 slice_axis=recon3d.types.CartesianAxis3D.Z,
922 parent_dir=save_path,
923 folder_name=instance_volumes.name,
924 file_type=".tif",
925 pad_length=4,
926 )
929def semantic_to_instance(
930 semantic_stack: SemanticImageStack,
931 instance_name: str,
932 instance_value: int,
933 min_feature_size: int,
934) -> InstanceImageStack:
935 """
936 Isolate each unique class in the semantic segmentation (e.g., Cats or Dogs)
937 and create the Instance Image Stacks.
939 Parameters
940 ----------
941 semantic_stack : SemanticImageStack
942 The semantic image stack containing the segmented data.
943 instance_name : str
944 The name for the instance image stack.
945 instance_value : int
946 The value in the semantic stack representing the class to isolate.
947 min_feature_size : int
948 The minimum feature size to retain in the instance image stack.
950 Returns
951 -------
952 InstanceImageStack
953 The instance image stack with isolated features.
955 Examples
956 --------
957 >>> semantic_stack = SemanticImageStack(
958 ... name="example_semantic_stack",
959 ... metadata=MetaData(
960 ... data_volume=DataVolume(z_image_count=10, y_height=256, x_width=256),
961 ... resolution=Resolution(
962 ... dx=Length(1.0, Units.MICRON),
963 ... dy=Length(1.0, Units.MICRON),
964 ... dz=Length(1.0, Units.MICRON)
965 ... ),
966 ... pixel_units=Units.MICRON,
967 ... origin=Origin(
968 ... x0=Length(0.0, Units.MICRON),
969 ... y0=Length(0.0, Units.MICRON),
970 ... z0=Length(0.0, Units.MICRON)
971 ... )
972 ... ),
973 ... data=np.random.randint(0, 2, (10, 256, 256))
974 ... )
975 >>> instance_stack = semantic_to_instance(
976 ... semantic_stack,
977 ... instance_name="example_instance_stack",
978 ... instance_value=1,
979 ... min_feature_size=10
980 ... )
981 >>> print(instance_stack)
982 InstanceImageStack(name='example_instance_stack', metadata=MetaData(data_volume=DataVolume(z_image_count=10, y_height=256, x_width=256), resolution=Resolution(dx=Length(value=1.0, unit=<Units.MICRON: 'micron'>), dy=Length(value=1.0, unit=<Units.MICRON: 'micron'>), dz=Length(value=1.0, unit=<Units.MICRON: 'micron'>)), pixel_units=<Units.MICRON: 'micron'>, origin=Origin(x0=Length(value=0.0, unit=<Units.MICRON: 'micron'>), y0=Length(value=0.0, unit=<Units.MICRON: 'micron'>), z0=Length(value=0.0, unit=<Units.MICRON: 'micron'>))), data=array(...), nlabels=..., min_feature_size=10)
983 """
985 # Create a boolean mask of the array
986 masked_stack = semantic_stack.data == instance_value
988 print(f"\tlabelling connected components in '{instance_name}'")
990 # Use cc3d for connected components labeling
991 cc3d_instance_stack, cc3d_nlabels = cc3d.connected_components(
992 masked_stack,
993 connectivity=cs.CC3DConnectivity.FACE_NEIGHBORS,
994 return_N=True,
995 )
996 print(
997 f"\t\twith cc3d package, found {cc3d_nlabels} connected components in '{instance_name}'"
998 )
999 return InstanceImageStack(
1000 name=instance_name,
1001 data=cc3d_instance_stack,
1002 metadata=semantic_stack.metadata,
1003 nlabels=cc3d_nlabels,
1004 min_feature_size=min_feature_size,
1005 )
1008##### MAIN FUNCTIONS
1011# TODO remove hardcoding here with yml config type
1012def instance_analysis_included(settings: dict, key: str) -> bool:
1013 """
1014 Checks YAML settings to determine if instance analysis should be performed.
1016 This function checks the provided settings dictionary to see if instance analysis
1017 is included for the specified key.
1019 Parameters
1020 ----------
1021 settings : dict
1022 The settings dictionary, typically loaded from a YAML file.
1023 key : str
1024 The key representing the class label to check for instance analysis inclusion.
1026 Returns
1027 -------
1028 bool
1029 True if instance analysis is included for the specified key, False otherwise.
1031 Examples
1032 --------
1033 >>> settings = {
1034 ... "class_labels": {
1035 ... "cat": {
1036 ... "instance_analysis": {
1037 ... "include": True
1038 ... }
1039 ... },
1040 ... "dog": {
1041 ... "instance_analysis": {
1042 ... "include": False
1043 ... }
1044 ... }
1045 ... }
1046 ... }
1047 >>> instance_analysis_included(settings, "cat")
1048 True
1049 >>> instance_analysis_included(settings, "dog")
1050 False
1051 >>> instance_analysis_included(settings, "bird")
1052 False
1053 """
1055 labels = settings["class_labels"]
1056 try:
1057 return labels[key]["instance_analysis"]["include"]
1058 except KeyError:
1059 return False
1062def instance_properties(
1063 instance_stack: InstanceImageStack,
1064 inst_indices: InstanceIndices,
1065) -> InstanceProperties:
1066 """
1067 Calculate various properties for each instance in an InstanceImageStack.
1069 Parameters
1070 ----------
1071 instance_stack : InstanceImageStack
1072 The instance image stack containing the segmented data.
1073 inst_indices : InstanceIndices
1074 The instance indices containing the indices of each instance.
1076 Returns
1077 -------
1078 InstanceProperties
1079 The properties of each instance, including number of voxels, equivalent spherical diameters,
1080 centroids, ellipsoids, surface areas, and volumes.
1082 Examples
1083 --------
1084 >>> instance_stack = InstanceImageStack(
1085 ... name="example_instance_stack",
1086 ... metadata=MetaData(
1087 ... data_volume=DataVolume(z_image_count=10, y_height=256, x_width=256),
1088 ... resolution=Resolution(
1089 ... dx=Length(1.0, Units.MICRON),
1090 ... dy=Length(1.0, Units.MICRON),
1091 ... dz=Length(1.0, Units.MICRON)
1092 ... ),
1093 ... pixel_units=Units.MICRON,
1094 ... origin=Origin(
1095 ... x0=Length(0.0, Units.MICRON),
1096 ... y0=Length(0.0, Units.MICRON),
1097 ... z0=Length(0.0, Units.MICRON)
1098 ... )
1099 ... ),
1100 ... data=np.random.randint(0, 2, (10, 256, 256)),
1101 ... nlabels=3,
1102 ... min_feature_size=2
1103 ... )
1104 >>> inst_indices = InstanceIndices(
1105 ... source_name="example_instance_stack",
1106 ... indices=[np.array([[0, 0], [1, 1]]), np.array([[0, 1], [1, 2]])],
1107 ... labels=InstanceLabels(data=np.array([1, 2]))
1108 ... )
1109 >>> properties = instance_properties(instance_stack, inst_indices)
1110 >>> print(properties)
1111 InstanceProperties(source_name='example_instance_stack', labels=InstanceLabels(data=array([1, 2])), n_voxels=[NVoxel(value=0), NVoxel(value=2), NVoxel(value=2)], equivalent_sphere_diameters=[Length(value=0.0, unit=<Units.MICRON: 'micron'>), Length(value=1.240700981798799, unit=<Units.MICRON: 'micron'>), Length(value=1.240700981798799, unit=<Units.MICRON: 'micron'>)], centroids=Centroids(data=[Centroid(cx=Length(value=0.0, unit=<Units.MICRON: 'micron'>), cy=Length(value=0.0, unit=<Units.MICRON: 'micron'>), cz=Length(value=0.0, unit=<Units.MICRON: 'micron'>)), Centroid(cx=Length(value=0.5, unit=<Units.MICRON: 'micron'>), cy=Length(value=0.5, unit=<Units.MICRON: 'micron'>), cz=Length(value=0.5, unit=<Units.MICRON: 'micron'>)), Centroid(cx=Length(value=1.5, unit=<Units.MICRON: 'micron'>), cy=Length(value=1.5, unit=<Units.MICRON: 'micron'>), cz=Length(value=1.5, unit=<Units.MICRON: 'micron'>))]), ellipsoids=BestFitEllipsoids(data=[BestFitEllipsoid(a=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=0.0, w=0.0)), b=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=0.0, w=0.0)), c=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=0.0, w=0.0))]), BestFitEllipsoid(a=EllipsoidAxis(length=Length(value=1.240700981798799, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=1.0, v=0.0, w=0.0)), b=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=1.0, w=0.0)), c=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=0.0, w=1.0))]), BestFitEllipsoid(a=EllipsoidAxis(length=Length(value=1.240700981798799, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=1.0, v=0.0, w=0.0)), b=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=1.0, w=0.0)), c=EllipsoidAxis(length=Length(value=0.0, unit=<Units.MICRON: 'micron'>), orientation=UnitVector(u=0.0, v=0.0, w=1.0))])]), surface_areas=EllipsoidSurfaceAreas(data=[Area(value=0.0, unit_squared=<Units.MICRON: 'micron'>), Area(value=19.739208802178716, unit_squared=<Units.MICRON: 'micron'>), Area(value=19.739208802178716, unit_squared=<Units.MICRON: 'micron'>)]), volumes=EllipsoidVolumes(data=[Volume(value=0.0, unit_cubed=<Units.MICRON: 'micron'>), Volume(value=4.1887902047863905, unit_cubed=<Units.MICRON: 'micron'>), Volume(value=4.1887902047863905, unit_cubed=<Units.MICRON: 'micron'>)]))
1112 """
1114 num_labels = instance_stack.nlabels
1115 # feature 0 is the background, and we want to include the last label
1116 stack_n_voxels = num_voxels(indices=inst_indices)
1117 stack_eq_diameters = np.zeros(num_labels + 1, dtype=object)
1118 stack_centroids = np.zeros(num_labels + 1, dtype=object)
1119 stack_ellipsoids = np.zeros(num_labels + 1, dtype=object)
1120 stack_surface_areas = np.zeros(num_labels + 1, dtype=object)
1121 stack_volumes = np.zeros(num_labels + 1, dtype=object)
1123 resolution = instance_stack.metadata.resolution
1124 origin = instance_stack.metadata.origin
1126 # null types
1127 null_length = Length(value=0.0, unit=resolution.dx.unit)
1128 null_area = Area(value=0.0, unit_squared=resolution.dx.unit)
1129 null_volume = Volume(value=0.0, unit_cubed=resolution.dx.unit)
1130 null_axis = EllipsoidAxis(
1131 length=null_length,
1132 orientation=UnitVector(
1133 u=0.0,
1134 v=0.0,
1135 w=0.0,
1136 ),
1137 )
1138 # set first entry to correct type of zero for feature 0
1139 stack_n_voxels[0] = NVoxel(value=0)
1140 stack_eq_diameters[0] = null_length
1141 stack_centroids[0] = Centroid(
1142 cx=null_length,
1143 cy=null_length,
1144 cz=null_length,
1145 )
1146 stack_ellipsoids[0] = BestFitEllipsoid(
1147 a=null_axis,
1148 b=null_axis,
1149 c=null_axis,
1150 )
1151 stack_surface_areas[0] = null_area
1152 stack_volumes[0] = null_volume
1154 for instance in range(1, num_labels + 1):
1155 if ((instance - 1) % 1000) == 0:
1156 print(f"\t\tAnalyzing feature {instance} of {num_labels}")
1157 indices = inst_indices.indices[instance]
1158 n_voxels = stack_n_voxels[instance].value
1160 # Equivalent spherical diameter
1161 stack_eq_diameters[instance] = equivalent_spherical_diameter(
1162 n_voxels=n_voxels,
1163 resolution=resolution,
1164 )
1166 # Centroid (real-space coordinates)
1167 centroid = center_of_mass(
1168 indices=indices,
1169 resolution=resolution,
1170 origin=origin,
1171 )
1172 stack_centroids[instance] = centroid
1174 # Ellipsoid properties
1175 ellipsoid = fit_ellipsoid(
1176 indices=indices,
1177 centroid=centroid,
1178 resolution=resolution,
1179 origin=origin,
1180 )
1181 stack_ellipsoids[instance] = ellipsoid
1183 surface_area = ellipsoid_surface_area(ellipsoid=ellipsoid)
1184 stack_surface_areas[instance] = surface_area
1186 volume = ellipsoid_volume(ellipsoid=ellipsoid)
1187 stack_volumes[instance] = volume
1189 properties = InstanceProperties(
1190 source_name=instance_stack.name,
1191 labels=inst_indices.labels,
1192 n_voxels=stack_n_voxels,
1193 equivalent_sphere_diameters=stack_eq_diameters.tolist(),
1194 centroids=Centroids(data=stack_centroids.tolist()),
1195 ellipsoids=BestFitEllipsoids(data=stack_ellipsoids.tolist()),
1196 surface_areas=EllipsoidSurfaceAreas(data=stack_surface_areas.tolist()),
1197 volumes=EllipsoidVolumes(data=stack_volumes.tolist()),
1198 )
1200 return properties
1203def process(yml_file: Path) -> bool:
1204 """
1205 Processes a directory of images specified in a YAML file, creates a semantic image stack,
1206 and performs instance analysis on all semantic label classes within the stack.
1207 Saves all instance stacks and properties into an HDF5 file containing the original semantic image stack.
1209 Parameters
1210 ----------
1211 yml_file : Path
1212 The fully pathed YAML file that specifies settings.
1214 Returns
1215 -------
1216 bool
1217 True if the process was successful, False otherwise.
1219 Examples
1220 --------
1221 >>> yml_file = Path("path/to/input.yml")
1222 >>> process(yml_file)
1223 True
1224 """
1226 semantic_stack = process_image_stack(yml_file)
1227 n_semantic_labels = np.unique(semantic_stack.data)
1229 # create h5
1230 h5_path = hio.image_to_voxel(yml_file)
1232 # start instance analysis
1233 db = ut.yaml_to_dict(yml_file)
1234 for instance_name in db["class_labels"]:
1236 # # check yml values match the labels in the semantic image stack
1237 # for semantic_label in n_semantic_labels:
1239 if not instance_analysis_included(settings=db, key=instance_name):
1240 continue
1242 print(f"Analyzing '{instance_name}' semantic label")
1244 # create the instance_stacks
1245 inst_stack = semantic_to_instance(
1246 semantic_stack=semantic_stack,
1247 instance_name=instance_name,
1248 instance_value=db["class_labels"][instance_name]["value"],
1249 min_feature_size=db["class_labels"][instance_name]["instance_analysis"][
1250 "min_feature_size"
1251 ],
1252 )
1254 # get instance_indices
1255 inst_indices = instance_indices(
1256 instance_stack=inst_stack,
1257 )
1258 if inst_stack.min_feature_size > 0:
1259 # overwrites instance_stack and inst_indices
1260 inst_stack, inst_indices = minimum_size_filter(
1261 initial_stack=inst_stack,
1262 initial_indices=inst_indices,
1263 )
1265 # add to h5
1266 hio.add_to_h5(
1267 inst_stack,
1268 h5_path=h5_path,
1269 h5_group="VoxelData",
1270 )
1271 print(f"\tlabelled '{instance_name}' voxel data added")
1273 # TODO not implemented instance indices yet
1274 # hio.add_to_h5(
1275 # inst_indices,
1276 # h5_path=h5_path,
1277 # h5_group=f"{instance_name}",
1278 # )
1280 print(f"\tcalculating '{instance_name}' instance properties...")
1281 # calculate indepedendent feature properties
1282 inst_properties = instance_properties(
1283 instance_stack=inst_stack,
1284 inst_indices=inst_indices,
1285 )
1287 # add to h5
1288 hio.add_to_h5(
1289 inst_properties,
1290 h5_path=h5_path,
1291 h5_group=f"{instance_name}",
1292 )
1293 print(f"\t'{instance_name}' instance property data added")
1295 nearest_neighbors = nearest_neighbor_distance(
1296 centroids=inst_properties.centroids
1297 )
1298 hio.add_to_h5(
1299 nearest_neighbors,
1300 h5_path=h5_path,
1301 h5_group=f"{instance_name}",
1302 )
1303 print(f"\t'{instance_name}' neighborhood property data added")
1305 return True
1308def main():
1309 """
1310 Command-line interface for running instance analysis on image stacks.
1312 This function serves as the entry point for terminal-based access to the instance analysis module.
1313 It is invoked from the command line using the 'instance_analysis' command specified in pyproject.toml.
1314 The function processes a YAML input file to create a semantic image stack and performs instance analysis
1315 on all semantic label classes within the stack. The results, including instance stacks and properties,
1316 are saved into an HDF5 file.
1318 Parameters
1319 ----------
1320 None
1322 Returns
1323 -------
1324 None
1326 Examples
1327 --------
1328 To run the instance analysis, use the following command in the terminal:
1329 $ instance_analysis path/to/input.yml
1330 """
1332 parser = argparse.ArgumentParser()
1333 parser.add_argument("input_file", help="the .yml user input file")
1334 args = parser.parse_args()
1335 input_file = args.input_file
1336 process(input_file)
1338 print(f"\n{input_file} processed!")
1341if __name__ == "__main__":
1342 main()