Coverage for src/recon3d/instance_analysis.py: 85%

1""" 

2Basic workflow for processing grayscale image stacks to labeled images for 

3feature property measurement. 

4 

5This module supports multiple phases (>2 phases) and includes simple 

6segmentation for demonstration purposes. 

7 

8It writes outputs as image directories, but can be integrated with HDF for 

9internal dataset storage. 

10""" 

11 

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 

18 

19import argparse 

20import glob 

21from typing import Final, NamedTuple, Tuple, Type, List 

22import math 

23 

24import h5py 

25import numpy as np 

26import numpy.typing as npt 

27from pathlib import Path 

28from PIL import Image 

29import skimage.io as skio 

30from sklearn.neighbors import NearestNeighbors as sklearnNN 

31 

32from scipy import ndimage 

33from skimage.measure import label, regionprops 

34import cc3d 

35import fastremap 

36 

37import recon3d.constants as cs 

38import recon3d.types 

39import recon3d.utility as ut 

40import recon3d.hdf_io as hio 

41from recon3d.types import * 

42 

43 

44# # ---------------------------------- 

45# # USER INPUT DATA HARD CODES - begin 

46# # ---------------------------------- 

47# # NEED TO PUT INTO A .yml file instead of hard codes 

48# # input_dir: str = "/Users/apolon/Desktop/test_recon/input/" 

49# # save_dir: str = "/Users/apolon/Desktop/test_recon/output/" 

50 

51# input_dir: Final[str] = "/Users/chovey/temp/test_recon/input/" 

52# save_dir: Final[str] = "/Users/chovey/temp/test_recon/output/" 

53 

54# file_type: Final[str] = ".tif" 

55# gray_path: Final[str] = f"{save_dir}/gray/" 

56# threshold: Final[int] = 185 

57# segmented_path: Final[str] = f"{save_dir}/segmented/" 

58# # first phase is always assumed to be the sample (1 or 255 in segmented images) 

59# phase_paths: np.ndarray = np.array( 

60# [f"{save_dir}/sample_mask/", f"{save_dir}/void_mask/"] 

61# ) 

62# # first label path is always assumed to the first phase after the sample, 

63# # same ordering as phase_paths 

64# label_paths: np.ndarray = np.array([f"{save_dir}/void_ids/"]) 

65 

66# # -------------------------------- 

67# # USER INPUT DATA HARD CODES - end 

68# # -------------------------------- 

69 

70 

71def calc_moment( 

72 indices: np.ndarray[np.int32], 

73 p: int, 

74 q: int, 

75 r: int, 

76 centroid_px: Centroid, 

77) -> float: 

78 """ 

79 Calculate 3D moment invariants of a feature based on its voxel locations. 

80 Mask is given in index locations, so math here is all done in pixel units. 

81 

82 Parameters 

83 ---------- 

84 indices : np.ndarray[np.int32] 

85 Array of [Z,Y,X] indices corresponding to the feature's voxel 

86 locations. 

87 The shape of the array should be (n, 3), where n is the number of 

88 voxels of the feature. 

89 Example: np.array([[Z1, Y1, X1], [Z2, Y2, X2], ..., [Zn, Yn, Xn]]). 

90 p : int 

91 Order of the moment for the X dimension. 

92 q : int 

93 Order of the moment for the Y dimension. 

94 r : int 

95 Order of the moment for the Z dimension. 

96 centroid_px : Centroid 

97 Centroid location for X, Y, and Z axis, respectively, in pixel units. 

98 

99 Returns 

100 ------- 

101 float 

102 The 3D central moment (translation invariant) in pixel units 

103 normalized by object volume. 

104 

105 Raises 

106 ------ 

107 ValueError 

108 If the units of the centroid are not in pixel/voxel units. 

109 

110 Examples 

111 -------- 

112 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32) 

113 >>> centroid_px = Centroid( 

114 ... cx=Length(1.0, Units.VOXEL), 

115 ... cy=Length(1.0, Units.VOXEL), 

116 ... cz=Length(1.0, Units.VOXEL) 

117 ... ) 

118 >>> calc_moment(indices, 2, 2, 2, centroid_px) 

119 0.66 

120 """ 

121 

122 unit = centroid_px.cx.unit 

123 if unit != Units.VOXEL: 

124 raise ValueError( 

125 f"Moment invariants are calculated using pixel/voxel units, but {unit} units were given" 

126 ) 

127 

128 mu = 0 

129 for i in range(0, indices.shape[0]): 

130 mu += ( 

131 (indices[i][2] - centroid_px.cx.value) 

132 ** p # x-location of voxel w.r.t x-centroid 

133 * (indices[i][1] - centroid_px.cy.value) 

134 ** q # y-location of voxel w.r.t y-centroid 

135 * (indices[i][0] - centroid_px.cz.value) 

136 ** r # z-location of voxel w.r.t z-centroid 

137 ) 

138 return mu / len(indices) 

139 

140 

141def center_of_mass( 

142 indices: np.ndarray[np.int32], resolution: Resolution, origin: Origin 

143) -> Centroid: 

144 """ 

145 Calculates the center of mass (centroid) of an instance object. 

146 Requires that the instance object has uniform density (is homogeneous). 

147 

148 Parameters 

149 ---------- 

150 indices : np.ndarray[np.int32] 

151 Array of [Z,Y,X] indices corresponding to the feature's voxel 

152 locations. 

153 The shape of the array should be (n, 3), where n is the number of 

154 voxels of the feature. 

155 Example: np.array([[Z1, Y1, X1], [Z2, Y2, X2], ..., [Zn, Yn, Xn]]). 

156 resolution : Resolution 

157 Resolution of the voxel grid in each dimension (dx, dy, dz). 

158 origin : Origin 

159 Origin of the voxel grid in each dimension (x0, y0, z0). 

160 

161 Returns 

162 ------- 

163 Centroid 

164 The centroid location for X, Y, and Z axis, respectively, in microns. 

165 

166 Examples 

167 -------- 

168 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32) 

169 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON)) 

170 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON)) 

171 >>> center_of_mass(indices, resolution, origin) 

172 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)) 

173 """ 

174 

175 centroid_zyx = np.add( 

176 np.multiply( 

177 np.mean(indices, axis=0), 

178 np.array( 

179 [ 

180 resolution.dz.value, 

181 resolution.dy.value, 

182 resolution.dx.value, 

183 ] 

184 ), 

185 ), 

186 np.array( 

187 [ 

188 origin.z0.value, 

189 origin.y0.value, 

190 origin.x0.value, 

191 ] 

192 ), 

193 ) 

194 cx = centroid_zyx[2] 

195 cy = centroid_zyx[1] 

196 cz = centroid_zyx[0] 

197 unit = resolution.dx.unit 

198 return Centroid( 

199 cx=Length(cx, unit=unit), 

200 cy=Length(cy, unit=unit), 

201 cz=Length(cz, unit=unit), 

202 ) 

203 

204 

205def centroid_pix( 

206 centroid: Centroid, resolution: Resolution, origin: Origin 

207) -> Centroid: 

208 """ 

209 Convert between real space units and pixel units. 

210 

211 Parameters 

212 ---------- 

213 centroid : Centroid 

214 The centroid location in real space units. 

215 resolution : Resolution 

216 Resolution of the voxel grid in each dimension (dx, dy, dz). 

217 origin : Origin 

218 Origin of the voxel grid in each dimension (x0, y0, z0). 

219 

220 Returns 

221 ------- 

222 Centroid 

223 The centroid location in pixel units. 

224 

225 Examples 

226 -------- 

227 >>> centroid = Centroid(cx=Length(10.0, Units.MICRON), cy=Length(20.0, Units.MICRON), cz=Length(30.0, Units.MICRON)) 

228 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON)) 

229 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON)) 

230 >>> centroid_pix(centroid, resolution, origin) 

231 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)) 

232 """ 

233 

234 # convert centroids to normalized pixel units 

235 cx = (centroid.cx.value - origin.x0.value) / resolution.dx.value 

236 cy = (centroid.cy.value - origin.y0.value) / resolution.dy.value 

237 cz = (centroid.cz.value - origin.z0.value) / resolution.dz.value 

238 

239 unit = Units.VOXEL 

240 

241 return Centroid( 

242 cx=Length(cx, unit=unit), 

243 cy=Length(cy, unit=unit), 

244 cz=Length(cz, unit=unit), 

245 ) 

246 

247 

248def ellipsoid_surface_area(ellipsoid: BestFitEllipsoid) -> Area: 

249 """ 

250 Calculate the surface area of the best fit ellipsoid using the Knud Thomsen formula. 

251 Applicable to scalene ellipsoids (semi-axes a > b > c), with a relative error of at most 1.061%. 

252 

253 Parameters 

254 ---------- 

255 ellipsoid : BestFitEllipsoid 

256 The best fit ellipsoid with semi-axes lengths a, b, and c. 

257 

258 Returns 

259 ------- 

260 Area 

261 The surface area of the ellipsoid in microns squared. 

262 

263 Examples 

264 -------- 

265 >>> ellipsoid = BestFitEllipsoid( 

266 ... a=EllipsoidAxis(length=Length(5.0, Units.MICRON), orientation=UnitVector(u=1, v=0, w=0)), 

267 ... b=EllipsoidAxis(length=Length(3.0, Units.MICRON), orientation=UnitVector(u=0, v=1, w=0)), 

268 ... c=EllipsoidAxis(length=Length(2.0, Units.MICRON), orientation=UnitVector(u=0, v=0, w=1)) 

269 ... ) 

270 >>> ellipsoid_surface_area(ellipsoid) 

271 Area(value=134.8149867941625, unit_squared=<Units.MICRON: 'micron'>) 

272 """ 

273 unit = ellipsoid.a.length.unit 

274 p = cs.Constants().KNUD_THOMSEN_FACTOR 

275 a, b, c = ( 

276 ellipsoid.a.length.value, 

277 ellipsoid.b.length.value, 

278 ellipsoid.c.length.value, 

279 ) 

280 

281 surface_area = ( 

282 4 

283 * math.pi 

284 * ((((a**p * b**p) + (a**p * c**p) + (b**p * c**p)) / 3) ** (1 / p)) 

285 ) 

286 return Area(value=surface_area, unit_squared=unit) 

287 

288 

289def ellipsoid_volume(ellipsoid: BestFitEllipsoid) -> Volume: 

290 """ 

291 Calculate the volume of the best fit ellipsoid using semi-axes lengths. 

292 

293 Parameters 

294 ---------- 

295 ellipsoid : BestFitEllipsoid 

296 The best fit ellipsoid with semi-axes lengths a, b, and c. 

297 

298 Returns 

299 ------- 

300 Volume 

301 The volume of the ellipsoid in microns cubed. 

302 

303 Examples 

304 -------- 

305 >>> ellipsoid = BestFitEllipsoid( 

306 ... a=EllipsoidAxis(length=Length(5.0, Units.MICRON), orientation=UnitVector(u=1, v=0, w=0)), 

307 ... b=EllipsoidAxis(length=Length(3.0, Units.MICRON), orientation=UnitVector(u=0, v=1, w=0)), 

308 ... c=EllipsoidAxis(length=Length(2.0, Units.MICRON), orientation=UnitVector(u=0, v=0, w=1)) 

309 ... ) 

310 >>> ellipsoid_volume(ellipsoid) 

311 Volume(value=125.66370614359171, unit_cubed=<Units.MICRON: 'micron'>) 

312 """ 

313 

314 unit = ellipsoid.a.length.unit 

315 a, b, c = ( 

316 ellipsoid.a.length.value, 

317 ellipsoid.b.length.value, 

318 ellipsoid.c.length.value, 

319 ) 

320 

321 volume = 4 / 3 * math.pi * a * b * c 

322 

323 return Volume(value=volume, unit_cubed=unit) 

324 

325 

326def equivalent_spherical_diameter( 

327 n_voxels: int, 

328 resolution: Resolution, 

329) -> Length: 

330 """ 

331 Calculate the equivalent spherical diameter of a feature. 

332 

333 Parameters 

334 ---------- 

335 n_voxels : int 

336 The number of voxels in the feature. 

337 resolution : Resolution 

338 The resolution of the voxel grid in each dimension (dx, dy, dz). 

339 

340 Returns 

341 ------- 

342 Length 

343 The equivalent spherical diameter of the feature in the same units as 

344 the resolution. 

345 

346 Examples 

347 -------- 

348 >>> resolution = Resolution( 

349 ... dx=Length(1.0, Units.MICRON), 

350 ... dy=Length(1.0, Units.MICRON), 

351 ... dz=Length(1.0, Units.MICRON) 

352 ... ) 

353 >>> equivalent_spherical_diameter(1000, resolution) 

354 Length(value=12.407009817988, unit=<Units.MICRON: 'micron'>) 

355 """ 

356 

357 value = 2 * np.power( 

358 ( 

359 (3 / (4 * np.pi)) 

360 * ( 

361 n_voxels 

362 * resolution.dx.value 

363 * resolution.dy.value 

364 * resolution.dz.value 

365 ) 

366 ), 

367 (1 / 3), 

368 ) 

369 return Length(value=value, unit=resolution.dx.unit) 

370 

371 

372def fit_ellipsoid( 

373 indices: np.ndarray, 

374 centroid: Centroid, 

375 resolution: Resolution, 

376 origin: Origin, 

377) -> BestFitEllipsoid: 

378 """ 

379 Fit a 3D ellipsoid to a set of voxel indices. 

380 

381 Parameters 

382 ---------- 

383 indices : np.ndarray 

384 Array of [Z,Y,X] indices corresponding to the feature's voxel 

385 locations. 

386 The shape of the array should be (n, 3), where n is the number of 

387 voxels of the feature. 

388 centroid : Centroid 

389 Centroid location for X, Y, and Z axis, respectively, in spatial 

390 coordinates (e.g., microns). 

391 resolution : Resolution 

392 Voxel resolution for X, Y, and Z dimensions, respectively, in length 

393 units (e.g., microns). 

394 origin : Origin 

395 Optional origin in spatial coordinates if the volume has been 

396 previously transformed. 

397 

398 Returns 

399 ------- 

400 BestFitEllipsoid 

401 The best fit ellipsoid with semi-axes lengths and orientations. 

402 

403 Examples 

404 -------- 

405 >>> indices = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]], dtype=np.int32) 

406 >>> centroid = Centroid(cx=Length(1.0, Units.MICRON), cy=Length(1.0, Units.MICRON), cz=Length(1.0, Units.MICRON)) 

407 >>> resolution = Resolution(dx=Length(1.0, Units.MICRON), dy=Length(1.0, Units.MICRON), dz=Length(1.0, Units.MICRON)) 

408 >>> origin = Origin(x0=Length(0.0, Units.MICRON), y0=Length(0.0, Units.MICRON), z0=Length(0.0, Units.MICRON)) 

409 >>> fit_ellipsoid(indices, centroid, resolution, origin) 

410 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))) 

411 """ 

412 

413 if not resolution.dx == resolution.dy == resolution.dz: 

414 raise NotImplementedError( 

415 "Only isotropic voxels are currently supported for axis length calculation" 

416 ) 

417 

418 unit = resolution.dx.unit 

419 

420 # convert centroids to normalized pixel units 

421 cen_px = centroid_pix(centroid, resolution, origin) 

422 

423 # Get moments 

424 s_xx = calc_moment(indices, 2, 0, 0, cen_px) 

425 s_yy = calc_moment(indices, 0, 2, 0, cen_px) 

426 s_zz = calc_moment(indices, 0, 0, 2, cen_px) 

427 s_xy = calc_moment(indices, 1, 1, 0, cen_px) 

428 s_xz = calc_moment(indices, 1, 0, 1, cen_px) 

429 s_yz = calc_moment(indices, 0, 1, 1, cen_px) 

430 

431 # Create matrix of 3D moments 

432 s_ij = np.array( 

433 [ 

434 [s_xx, s_xy, s_xz], 

435 [s_xy, s_yy, s_yz], 

436 [s_xz, s_yz, s_zz], 

437 ] 

438 ) 

439 eig_val, eig_vec = np.linalg.eig(s_ij) 

440 

441 semi_ax_lengths_px = np.sqrt(5 * eig_val) # still in pixel units 

442 sort = np.argsort(semi_ax_lengths_px)[::-1] # sort descending axis lengths 

443 

444 # TODO: math to account for change in length as function 

445 # as function of axis orientation in anisotropic voxels 

446 # remove NotImplementedError above 

447 semi_ax_lengths = np.array( 

448 [ 

449 semi_ax_lengths_px[sort[0]] * resolution.dx.value, 

450 semi_ax_lengths_px[sort[1]] * resolution.dx.value, 

451 semi_ax_lengths_px[sort[2]] * resolution.dx.value, 

452 ] 

453 ) 

454 # some axes lengths are NAN, should overwrite to 0.0 

455 semi_ax_lengths = np.nan_to_num(semi_ax_lengths, nan=0.0) 

456 

457 # from docs on np.linalg.eig: The normalized (unit “length”) eigenvectors, 

458 # such that the column eigenvectors[:,i] is the eigenvector corresponding to the eigenvalue eigenvalues[i]. 

459 ax_vectors = np.array( 

460 [ 

461 eig_vec[:, sort[0]], 

462 eig_vec[:, sort[1]], 

463 eig_vec[:, sort[2]], 

464 ] 

465 ) 

466 ax_vectors = np.reshape(ax_vectors, (9)) 

467 

468 a_vector = UnitVector( 

469 u=ax_vectors[0], 

470 v=ax_vectors[1], 

471 w=ax_vectors[2], 

472 ) 

473 a = EllipsoidAxis( 

474 length=Length(semi_ax_lengths[0], unit=unit), 

475 orientation=a_vector, 

476 ) 

477 b_vector = UnitVector( 

478 u=ax_vectors[3], 

479 v=ax_vectors[4], 

480 w=ax_vectors[5], 

481 ) 

482 b = EllipsoidAxis( 

483 length=Length(semi_ax_lengths[1], unit=unit), 

484 orientation=b_vector, 

485 ) 

486 c_vector = UnitVector( 

487 u=ax_vectors[6], 

488 v=ax_vectors[7], 

489 w=ax_vectors[8], 

490 ) 

491 c = EllipsoidAxis( 

492 length=Length(semi_ax_lengths[2], unit=unit), 

493 orientation=c_vector, 

494 ) 

495 

496 ellipsoid = BestFitEllipsoid(a=a, b=b, c=c) 

497 return ellipsoid 

498 

499 

500def instance_indices( 

501 instance_stack: InstanceImageStack, 

502) -> InstanceIndices: 

503 """ 

504 Returns an array of instances and a list of N-dimensional arrays of indices 

505 of the unique values in the instance stack. 

506 

507 Parameters 

508 ---------- 

509 instance_stack : InstanceImageStack 

510 Array with arbitrary dimensions. 

511 

512 Returns 

513 ------- 

514 InstanceIndices 

515 - 1D-array of sorted unique values. 

516 - List of arrays. Each array contains the indices where a given value in x is found. 

517 

518 Examples 

519 -------- 

520 >>> instance_stack = InstanceImageStack( 

521 ... name="example_stack", 

522 ... metadata=MetaData( 

523 ... data_volume=DataVolume(x_width=2, y_height=2, z_image_count=1), 

524 ... resolution=Resolution( 

525 ... dx=Length(1.0, Units.MICRON), 

526 ... dy=Length(1.0, Units.MICRON), 

527 ... dz=Length(1.0, Units.MICRON) 

528 ... ), 

529 ... pixel_units=Units.MICRON, 

530 ... origin=Origin( 

531 ... x0=Length(0.0, Units.MICRON), 

532 ... y0=Length(0.0, Units.MICRON), 

533 ... z0=Length(0.0, Units.MICRON) 

534 ... ) 

535 ... ), 

536 ... data=np.array([[1, 2, 2], [3, 1, 1]]), 

537 ... nlabels=3, 

538 ... min_feature_size=1 

539 ... ) 

540 >>> instance_indices(instance_stack) 

541 InstanceIndices( 

542 source_name='example_stack', 

543 indices=[array([[0, 0], [1, 1], [1, 2]]), array([[0, 1], [0, 2]]), array([[1, 0]])], 

544 labels=InstanceLabels(data=np.array([1, 2, 3])) 

545 ) 

546 """ 

547 

548 x = instance_stack.data 

549 x_flat = x.ravel() 

550 ix_flat = np.argsort(x_flat) 

551 u, ix_u = fastremap.unique(x_flat[ix_flat], return_index=True) 

552 ix_ndim = np.unravel_index(ix_flat, x.shape) 

553 ix_ndim = np.c_[ix_ndim] if x.ndim > 1 else ix_flat 

554 indices = np.split(ix_ndim, ix_u[1:]) 

555 labels = u 

556 instance_ids = InstanceIndices( 

557 source_name=instance_stack.name, 

558 indices=indices, 

559 labels=InstanceLabels(data=labels), 

560 ) 

561 return instance_ids 

562 

563 

564def minimum_size_filter( 

565 initial_stack: InstanceImageStack, 

566 initial_indices: InstanceIndices, 

567) -> Tuple[InstanceImageStack, InstanceIndices]: 

568 """ 

569 Removes features below minimum size from InstanceIndices and corresponding 

570 InstanceImageStack. Should be done before shape metrics are determined. 

571 

572 Parameters 

573 ---------- 

574 initial_stack : InstanceImageStack 

575 The initial instance image stack. 

576 initial_indices : InstanceIndices 

577 The initial instance indices. 

578 

579 Returns 

580 ------- 

581 Tuple[InstanceImageStack, InstanceIndices] 

582 The filtered instance image stack and the filtered instance indices. 

583 

584 Raises 

585 ------ 

586 ValueError 

587 If the names of the initial stack and indices do not match. 

588 

589 Examples 

590 -------- 

591 >>> initial_stack = InstanceImageStack( 

592 ... name="example_stack", 

593 ... metadata=MetaData( 

594 ... data_volume=DataVolume(x_width=4, y_height=2, z_image_count=1), 

595 ... resolution=Resolution( 

596 ... dx=Length(1.0, Units.MICRON), 

597 ... dy=Length(1.0, Units.MICRON), 

598 ... dz=Length(1.0, Units.MICRON) 

599 ... ), 

600 ... pixel_units=Units.MICRON, 

601 ... origin=Origin( 

602 ... x0=Length(0.0, Units.MICRON), 

603 ... y0=Length(0.0, Units.MICRON), 

604 ... z0=Length(0.0, Units.MICRON) 

605 ... ) 

606 ... ), 

607 ... data=np.array([[1, 2, 2, 1], [3, 1, 1, 1]]), 

608 ... nlabels=3, 

609 ... min_feature_size=2 

610 ... ) 

611 >>> initial_indices = instance_indices(initial_stack) 

612 >>> minimum_size_filter(initial_stack, initial_indices) 

613 2 connected components remaining after minimum size filter of 2 voxels 

614 (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]])])) 

615 """ 

616 

617 if not initial_stack.name == initial_indices.source_name: 

618 raise ValueError( 

619 f"""Initial InstanceImageStack name ({initial_stack.name}) and 

620 InstanceIndices name ({initial_indices.source_name}) do not match. 

621 Function inputs must be from same dataset.""" 

622 ) 

623 

624 n_voxels = num_voxels(initial_indices) 

625 n_voxels_vals = np.asarray([i.value for i in n_voxels]) 

626 invalid_instance_ids = np.nonzero(n_voxels_vals < initial_stack.min_feature_size)[0] 

627 masked_instance_ids = fastremap.mask(initial_stack.data, invalid_instance_ids) 

628 remapped_instance_ids, _ = fastremap.renumber( 

629 masked_instance_ids, 

630 preserve_zero=True, 

631 ) 

632 nlabels = np.max(np.unique(remapped_instance_ids)) 

633 

634 filtered_stack = InstanceImageStack( 

635 name=initial_stack.name, 

636 metadata=initial_stack.metadata, 

637 data=remapped_instance_ids, 

638 nlabels=nlabels, 

639 min_feature_size=initial_stack.min_feature_size, 

640 ) 

641 

642 filtered_indices = instance_indices(instance_stack=filtered_stack) 

643 

644 print( 

645 f"\t\t{nlabels} connected components remaining after minimum size filter of {initial_stack.min_feature_size} voxels" 

646 ) 

647 return filtered_stack, filtered_indices 

648 

649 

650# TODO 

651def map_features_to_voxels(): 

652 """To come.""" 

653 

654 

655def nearest_neighbor_distance( 

656 centroids: Centroids, 

657) -> NthNearestNeighbors: 

658 """ 

659 Calculate the nearest neighbor distance for a set of centroids. 

660 With nth_nearest_neighbor=1 returning the features themselves 

661 (distances=0, neighbor is themselves). 

662 

663 Parameters 

664 ---------- 

665 centroids : Centroids 

666 The centroids for which to calculate the nearest neighbor distance. 

667 

668 Returns 

669 ------- 

670 NthNearestNeighbors 

671 The nth nearest neighbors, distances, and instance IDs. 

672 

673 Examples 

674 -------- 

675 >>> centroids = Centroids(data=[ 

676 ... Centroid(cx=Length(0.0, Units.MICRON), cy=Length(0.0, Units.MICRON), cz=Length(0.0, Units.MICRON)), 

677 ... Centroid(cx=Length(1.0, Units.MICRON), cy=Length(1.0, Units.MICRON), cz=Length(1.0, Units.MICRON)), 

678 ... Centroid(cx=Length(2.0, Units.MICRON), cy=Length(2.0, Units.MICRON), cz=Length(2.0, Units.MICRON)) 

679 ... ]) 

680 >>> nearest_neighbor_distance(centroids) 

681 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])) 

682 """ 

683 

684 neighbor = sklearnNN(n_neighbors=2) 

685 centroid_array = ut.centroids_to_ndarray(centroids=centroids) 

686 neighbor.fit(centroid_array) 

687 neighbor_dist, neighbor_index = neighbor.kneighbors( 

688 centroid_array, 

689 return_distance=True, 

690 ) 

691 distances = neighbor_dist[:, 1] 

692 instance_id = neighbor_index[:, 1] 

693 distances[0], instance_id[0] = 0, 0 # set feature 0 to 0 distance and neighbor 0 

694 

695 unit = centroids.data[0].cx.unit 

696 distances_list = [Length(value=i, unit=unit) for i in distances.tolist()] 

697 nth_nearest = NthNearestNeighbors( 

698 nth_nearest=2, 

699 distances=distances_list, 

700 instance_id=instance_id, 

701 ) 

702 return nth_nearest 

703 

704 

705def num_voxels(indices: InstanceIndices) -> list[NVoxel]: 

706 """ 

707 Returns the number of voxels for each instance ID from an InstanceImageStack. 

708 

709 Parameters 

710 ---------- 

711 indices : InstanceIndices 

712 The instance indices containing the indices of each instance. 

713 

714 Returns 

715 ------- 

716 List[NVoxel] 

717 A list of NVoxel objects, each representing the number of voxels for an instance. 

718 

719 Examples 

720 -------- 

721 >>> indices = InstanceIndices( 

722 ... source_name="example_stack", 

723 ... indices=[np.array([[0, 0], [1, 1]]), np.array([[0, 1], [1, 2]])], 

724 ... labels=InstanceLabels(data=np.array([1, 2])) 

725 ... ) 

726 >>> num_voxels(indices) 

727 [NVoxel(value=2, unit=<Units.VOXEL: 'voxel'>), NVoxel(value=2, unit=<Units.VOXEL: 'voxel'>)] 

728 """ 

729 

730 indices = indices.indices 

731 data = [NVoxel(value=len(i)) for i in indices] 

732 # n_voxels = np.asarray(data, dtype=NVoxel) 

733 n_voxels = data 

734 return n_voxels 

735 

736 

737def process_image_stack(yml_input_file: Path) -> SemanticImageStack: 

738 """ 

739 Given a yml input file, converts it to a SemanticImageStack. 

740 

741 Parameters 

742 ---------- 

743 yml_input_file : Path 

744 The fully pathed yml file that specifies settings. 

745 

746 Returns 

747 ------- 

748 SemanticImageStack 

749 The SemanticImageStack object. 

750 

751 Examples 

752 -------- 

753 >>> yml_input_file = Path("path/to/input.yml") 

754 >>> semantic_image_stack = process_image_stack(yml_input_file) 

755 """ 

756 

757 print(f"Processing specification file: {yml_input_file}") 

758 db = ut.yaml_to_dict(yml_input_file) 

759 

760 file_type = db["semantic_images_type"] 

761 valid_extensions = (".tif", ".tiff") 

762 

763 if file_type not in valid_extensions: 

764 raise ValueError("Error: 'image_type' must be '.tif' or '.tiff'.") 

765 

766 path_input = Path(db["semantic_images_dir"]).expanduser() 

767 if not path_input.is_dir(): 

768 raise ValueError(f"Error, 'image_dir' of {path_input} not found.") 

769 print(f"Input path: {path_input}") 

770 

771 path_output = Path(db["out_dir"]).expanduser() 

772 if not path_output.is_dir(): 

773 Path(path_output).mkdir( 

774 parents=True, exist_ok=True 

775 ) # makes directory if it doesn't exist 

776 print(f"Output path: {path_output}") 

777 

778 # semantic_seg_dict = db["class_labels"] 

779 

780 # # load images 

781 # imageStack = skio.imread_collection(tiff_files, conserve_memory=True) 

782 data = ut.read_images(path_input, file_type) 

783 semantic_seg_stack_name = db["semantic_imagestack_name"] 

784 

785 # create meta data 

786 (nz, ny, nx) = data.shape 

787 data_volume = DataVolume(z_image_count=nz, y_height=ny, x_width=nx) 

788 print( 

789 f"image stack, {semantic_seg_stack_name}, has dimensions (num_images, row, col): {data.shape}" 

790 ) 

791 

792 pixel_units = db["pixel_units"] 

793 valid_units = set(item.value for item in Units) 

794 if pixel_units not in valid_units: 

795 raise ValueError( 

796 f"Error, '{pixel_units}' is not a valid unit, accepted units are: {valid_units}" 

797 ) 

798 pixel_units = Units(db["pixel_units"]) 

799 

800 dx, dy, dz = ( 

801 db["voxel_size"]["dx"], 

802 db["voxel_size"]["dy"], 

803 db["voxel_size"]["dz"], 

804 ) 

805 resolution = Resolution( 

806 dx=Length(dx, unit=pixel_units), 

807 dy=Length(dy, unit=pixel_units), 

808 dz=Length(dz, unit=pixel_units), 

809 ) 

810 

811 x0, y0, z0 = ( 

812 db["origin"]["x0"], 

813 db["origin"]["y0"], 

814 db["origin"]["z0"], 

815 ) 

816 origin = Origin( 

817 x0=Length(x0, unit=pixel_units), 

818 y0=Length(y0, unit=pixel_units), 

819 z0=Length(z0, unit=pixel_units), 

820 ) 

821 

822 meta = MetaData( 

823 data_volume=data_volume, 

824 resolution=resolution, 

825 pixel_units=pixel_units, 

826 origin=origin, 

827 ) 

828 

829 # # Optional, grey image stack 

830 # # TODO: test functionality 

831 # if db["grey_image_dir"]: 

832 # path_grey_image_dir = Path(db["path_grey_image_dir"]).expanduser() 

833 # assert path_grey_image_dir.is_dir(), "Error, 'image_dir' not found." 

834 # grey_image_name = "Raw greyscale" 

835 # grey_data = ut.read_images(path_grey_image_dir, file_type) 

836 # assert data.shape == grey_data.shape 

837 # semantic_seg_image_stack = ImageStack( 

838 # name=grey_image_name, metadata=meta, data=grey_data 

839 # ) 

840 

841 # create SemanticImageStack 

842 semantic_seg_image_stack = SemanticImageStack( 

843 name=semantic_seg_stack_name, metadata=meta, data=data 

844 ) 

845 

846 return semantic_seg_image_stack # , path_file_output) 

847 

848 

849def save_instance_images( 

850 input_volumes: list[InstanceImageStack], save_path: Path 

851) -> bool: 

852 """ 

853 Save InstanceImageStack objects to image files. 

854 

855 This function takes a list of InstanceImageStack objects and saves each one as a series of image files. 

856 The images are saved in the specified directory, with each stack saved in its own subdirectory. 

857 

858 Parameters 

859 ---------- 

860 input_volumes : list of InstanceImageStack 

861 List of InstanceImageStack objects to be saved as images. 

862 save_path : Path 

863 The directory where the images will be saved. 

864 

865 Returns 

866 ------- 

867 bool 

868 True if the images were saved successfully, False otherwise. 

869 

870 Examples 

871 -------- 

872 >>> input_volumes = [ 

873 ... InstanceImageStack( 

874 ... name="example_stack_1", 

875 ... metadata=MetaData( 

876 ... data_volume=DataVolume(x_width=256, y_height=256, z_image_count=10), 

877 ... resolution=Resolution( 

878 ... dx=Length(1.0, Units.MICRON), 

879 ... dy=Length(1.0, Units.MICRON), 

880 ... dz=Length(1.0, Units.MICRON) 

881 ... ), 

882 ... pixel_units=Units.MICRON, 

883 ... origin=Origin( 

884 ... x0=Length(0.0, Units.MICRON), 

885 ... y0=Length(0.0, Units.MICRON), 

886 ... z0=Length(0.0, Units.MICRON) 

887 ... ) 

888 ... ), 

889 ... data=np.random.rand(10, 256, 256), 

890 ... nlabels=3, 

891 ... min_feature_size=2 

892 ... ), 

893 ... InstanceImageStack( 

894 ... name="example_stack_2", 

895 ... metadata=MetaData( 

896 ... data_volume=DataVolume(x_width=256, y_height=256, z_image_count=10), 

897 ... resolution=Resolution( 

898 ... dx=Length(1.0, Units.MICRON), 

899 ... dy=Length(1.0, Units.MICRON), 

900 ... dz=Length(1.0, Units.MICRON) 

901 ... ), 

902 ... pixel_units=Units.MICRON, 

903 ... origin=Origin( 

904 ... x0=Length(0.0, Units.MICRON), 

905 ... y0=Length(0.0, Units.MICRON), 

906 ... z0=Length(0.0, Units.MICRON) 

907 ... ) 

908 ... ), 

909 ... data=np.random.rand(10, 256, 256), 

910 ... nlabels=3, 

911 ... min_feature_size=2 

912 ... ) 

913 ... ] 

914 >>> save_path = Path("path/to/save") 

915 >>> save_instance_images(input_volumes, save_path) 

916 True 

917 """ 

918 

919 for instance_volumes in input_volumes: 

920 ut.ndarray_to_img( 

921 data=instance_volumes.data, 

922 slice_axis=recon3d.types.CartesianAxis3D.Z, 

923 parent_dir=save_path, 

924 folder_name=instance_volumes.name, 

925 file_type=".tif", 

926 pad_length=4, 

927 ) 

928 

929 

930def semantic_to_instance( 

931 semantic_stack: SemanticImageStack, 

932 instance_name: str, 

933 instance_value: int, 

934 min_feature_size: int, 

935) -> InstanceImageStack: 

936 """ 

937 Isolate each unique class in the semantic segmentation (e.g., Cats or Dogs) 

938 and create the Instance Image Stacks. 

939 

940 Parameters 

941 ---------- 

942 semantic_stack : SemanticImageStack 

943 The semantic image stack containing the segmented data. 

944 instance_name : str 

945 The name for the instance image stack. 

946 instance_value : int 

947 The value in the semantic stack representing the class to isolate. 

948 min_feature_size : int 

949 The minimum feature size to retain in the instance image stack. 

950 

951 Returns 

952 ------- 

953 InstanceImageStack 

954 The instance image stack with isolated features. 

955 

956 Examples 

957 -------- 

958 >>> semantic_stack = SemanticImageStack( 

959 ... name="example_semantic_stack", 

960 ... metadata=MetaData( 

961 ... data_volume=DataVolume(z_image_count=10, y_height=256, x_width=256), 

962 ... resolution=Resolution( 

963 ... dx=Length(1.0, Units.MICRON), 

964 ... dy=Length(1.0, Units.MICRON), 

965 ... dz=Length(1.0, Units.MICRON) 

966 ... ), 

967 ... pixel_units=Units.MICRON, 

968 ... origin=Origin( 

969 ... x0=Length(0.0, Units.MICRON), 

970 ... y0=Length(0.0, Units.MICRON), 

971 ... z0=Length(0.0, Units.MICRON) 

972 ... ) 

973 ... ), 

974 ... data=np.random.randint(0, 2, (10, 256, 256)) 

975 ... ) 

976 >>> instance_stack = semantic_to_instance( 

977 ... semantic_stack, 

978 ... instance_name="example_instance_stack", 

979 ... instance_value=1, 

980 ... min_feature_size=10 

981 ... ) 

982 >>> print(instance_stack) 

983 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) 

984 """ 

985 

986 # Create a boolean mask of the array 

987 masked_stack = semantic_stack.data == instance_value 

988 

989 print(f"\tlabelling connected components in '{instance_name}'") 

990 

991 # Use cc3d for connected components labeling 

992 cc3d_instance_stack, cc3d_nlabels = cc3d.connected_components( 

993 masked_stack, 

994 connectivity=cs.CC3DConnectivity.FACE_NEIGHBORS, 

995 return_N=True, 

996 ) 

997 print( 

998 f"\t\twith cc3d package, found {cc3d_nlabels} connected components in '{instance_name}'" 

999 ) 

1000 return InstanceImageStack( 

1001 name=instance_name, 

1002 data=cc3d_instance_stack, 

1003 metadata=semantic_stack.metadata, 

1004 nlabels=cc3d_nlabels, 

1005 min_feature_size=min_feature_size, 

1006 ) 

1007 

1008 

1009##### MAIN FUNCTIONS 

1010 

1011 

1012# TODO remove hardcoding here with yml config type 

1013def instance_analysis_included(settings: dict, key: str) -> bool: 

1014 """ 

1015 Checks YAML settings to determine if instance analysis should be performed. 

1016 

1017 This function checks the provided settings dictionary to see if instance analysis 

1018 is included for the specified key. 

1019 

1020 Parameters 

1021 ---------- 

1022 settings : dict 

1023 The settings dictionary, typically loaded from a YAML file. 

1024 key : str 

1025 The key representing the class label to check for instance analysis inclusion. 

1026 

1027 Returns 

1028 ------- 

1029 bool 

1030 True if instance analysis is included for the specified key, False otherwise. 

1031 

1032 Examples 

1033 -------- 

1034 >>> settings = { 

1035 ... "class_labels": { 

1036 ... "cat": {