Coverage for /usr/local/lib/python3.8/site-packages/spam/helpers/imageManipulation.py: 93%

1# Library of SPAM functions for manipulating images 

2# Copyright (C) 2020 SPAM Contributors 

3# 

4# This program is free software: you can redistribute it and/or modify it 

5# under the terms of the GNU General Public License as published by the Free 

6# Software Foundation, either version 3 of the License, or (at your option) 

7# any later version. 

8# 

9# This program is distributed in the hope that it will be useful, but WITHOUT 

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

12# more details. 

13# 

14# You should have received a copy of the GNU General Public License along with 

15# this program. If not, see <http://www.gnu.org/licenses/>. 

16 

17import numpy 

18import progressbar 

19import spam.label 

20import tifffile 

21 

22 

23def stackToArray(prefix, suffix=".tif", stack=range(10), digits="05d"): 

24 """ 

25 Convert of stack of 2D sequential tif images into a 3D array. 

26 

27 Parameters 

28 ---------- 

29 prefix : string 

30 The common name of the 2D images files before the sequential number 

31 

32 suffix : string, default='.tif' 

33 The common name and extension of the 2D images after the sequential number 

34 

35 stack : sequence, default=range(10) 

36 The numbers of the slices with no formating (with no leading zeros) 

37 

38 digits : string, default='05d' 

39 The format (number of digits) of the numbers (add leading zeros). 

40 

41 Returns 

42 ------- 

43 numpy array 

44 The 3D image 

45 """ 

46 

47 def _slice_name(pr, st, di, su): 

48 return f"{pr}{st:{di}}{su}" 

49 

50 # Step 1 If nBytes is not defined: we open the first slice just for the dimensions 

51 slice_name = _slice_name(prefix, stack[0], digits, suffix) 

52 slice_im = tifffile.imread(slice_name) 

53 

54 # Step 2 create empty array of good size and type 

55 ny, nx = slice_im.shape 

56 nz = len(stack) 

57 im = numpy.zeros([nz, ny, nx], dtype=slice_im.dtype) 

58 

59 # Step 3 stack all the slices 

60 for i, s in enumerate(stack): 

61 slice_name = _slice_name(prefix, s, digits, suffix) 

62 im[i, :, :] = tifffile.imread(slice_name) 

63 

64 return im 

65 

66 

67def crop(im, boxSize, boxOrigin=None): 

68 """ 

69 This function crops an image using slicing. 

70 

71 Parameters 

72 ---------- 

73 im: array 

74 The image to crop. 

75 

76 boxSize: int 

77 The size in voxels of the crop box (from boxOrigin). If a int, the size is the same for each axis. If a sequence, ``boxSize`` should contain one value for each axis. 

78 

79 boxOrigin : int, default=None 

80 The coordinates in voxels of the origin of the crop box. If a int, the coordinates are the same for each axis. 

81 If a tuple, ``boxOrigin`` should contain one value for each axis. If ``None`` the image is cropped from its centre. 

82 

83 Returns 

84 ------- 

85 array 

86 The cropped image. 

87 """ 

88 

89 # get center of the image 

90 imCentre = [int(s) / 2 for s in im.shape] 

91 # get sizes 

92 (sz, sy, sx) = (boxSize, boxSize, boxSize) if isinstance(boxSize, int) else boxSize 

93 

94 # get box origin 

95 if boxOrigin is not None: 

96 (cz, cy, cx) = (boxOrigin, boxOrigin, boxOrigin) if isinstance(boxOrigin, int) else boxOrigin 

97 else: 

98 (cz, cy, cx) = ( 

99 imCentre[0] - sz // 2, 

100 imCentre[1] - sy // 2, 

101 imCentre[2] - sx // 2, 

102 ) 

103 

104 # test sizes 

105 if cz + sz > im.shape[0] or cy + sy > im.shape[1] or cx + sx > im.shape[2]: 

106 print("spam.helpers.imageManipulation.crop: box bigger than image.") 

107 print("exit function.") 

108 return -1 

109 

110 return im[int(cz) : int(cz + sz), int(cy) : int(cy + sy), int(cx) : int(cx + sx)] 

111 

112 

113def rescale(im, scale=(0, 1)): 

114 """ 

115 This function **rescales** the values of an image according to a scale 

116 and save it to as 4 bytes floats (float32). 

117 

118 Parameters 

119 ---------- 

120 im: array 

121 The image to rescale 

122 

123 scale : (float, float), default=(0 1) 

124 The min and max of the rescaled image 

125 

126 Returns 

127 ------- 

128 array, float 

129 The rescaled image. 

130 

131 Examples 

132 -------- 

133 >>> im = numpy.random.randn( 100, 100, 100 ).astype( '<f4' ) 

134 produce float32 array of positive and negative numbers 

135 >>> imRescaled = rescale( im, scale=[-1, 1] ) 

136 produce float32 array of numbers between -1 and 1 

137 

138 """ 

139 

140 im_max = float(im.max()) 

141 im_min = float(im.min()) 

142 

143 return (min(scale) + (max(scale) - min(scale)) * ((im.astype("<f4") - im_min) / (im_max - im_min))).astype("<f4") 

144 

145 

146def rescaleToInteger(im, nBytes=1, scale=None): 

147 """ 

148 This function **rescales** a 4 bytes float image values to a unsigned integers of ``nBytes``. 

149 

150 Parameters 

151 ---------- 

152 im: float32 numpy array 

153 The image to rescale 

154 

155 nBytes : int, default=1 

156 The number of bytes of the unsigned interger output. 

157 Possible values are power of 2 

158 

159 .. code-block:: text 

160 

161 reminder 

162 1 byte = 8 bits -> ouput from 0 to 255 

163 2 bytes = 16 bits -> ouput from 0 to 65 535 

164 4 bytes = 32 bits -> ouput from 0 to 4 294 967 295 

165 

166 scale : (float, float), default=None 

167 If None, the maximum and minimum use for the rescaling is the maximum and the minimum of the image 

168 

169 Returns 

170 ------- 

171 numpy array, uint 

172 The rescaled image 

173 

174 Examples 

175 -------- 

176 >>> im = numpy.random.randn( 100, 100, 100 ).astype( '<f4' ) 

177 produce float32 array of positive and negative numbers 

178 >>> imRescaled = rescaleToInteger( im, nBytes=4 ) 

179 produce uint32 array of numbers between 0 and 4 294 967 295 

180 

181 """ 

182 

183 nBytes = int(nBytes) 

184 if (nBytes & (nBytes - 1)) != 0: 

185 raise ValueError(f"nBytes should be a power of 2 ({nBytes} given)") 

186 

187 # check if float32 given 

188 if im.dtype is not numpy.dtype("float32"): 

189 raise ValueError(f"image should be encode in float32 ({im.dtype} given)") 

190 

191 if scale is None: 

192 # if no scale is given: it takes the max and min of the image 

193 im_max = im.max() 

194 im_min = im.min() 

195 else: 

196 # if a scale is given take it if larger (smaller) than max (min) of image 

197 # im_max = max(scale) if max(scale) > im.max() else im.max() 

198 # im_min = min(scale) if min(scale) < im.min() else im.min() 

199 im_max = max(scale) 

200 im_min = min(scale) 

201 im[im > im_max] = im_max 

202 im[im < im_min] = im_min 

203 

204 im_min = float(im_min) 

205 im_max = float(im_max) 

206 

207 return ((2 ** (8 * nBytes) - 1) * ((im.astype("<f4") - im_min) / (im_max - im_min))).astype(f"<u{nBytes}") 

208 

209 

210def convertUnsignedIntegers(im, nBytes=1): 

211 """ 

212 This function **converts** an images of unsigned integers. 

213 

214 Note: this function does not rescale. 

215 

216 Parameters 

217 ---------- 

218 im: array, uint 

219 The image to convert. 

220 

221 nBytes : int, default=1 

222 The number of bytes of the unsigned interger output. 

223 Possible values are power of 2. 

224 

225 .. code-block:: text 

226 

227 reminder 

228 1 byte = 8 bits -> ouput from 0 to 255 

229 2 bytes = 16 bits -> ouput from 0 to 65 535 

230 4 bytes = 32 bits -> ouput from 0 to 4 294 967 295 

231 

232 Returns 

233 ------- 

234 array, uint 

235 The converted image. 

236 

237 Examples 

238 -------- 

239 >>> im = numpy.random.randint( 12, high=210, size=(100, 100, 100) ).astype( '<u1' ) 

240 produce an uint8 array of numbers between 12 and 210 

241 >>> imRescaled = rescaleToInteger( im, nBytes=2 ) 

242 produce an uint16 array 3084 and 53970 

243 

244 """ 

245 

246 nBytes = int(nBytes) 

247 if (nBytes & (nBytes - 1)) != 0: 

248 raise ValueError(f"nBytes should be a power of 2 ({nBytes} given)") 

249 

250 # number of bits of the output 

251 nbo = 8 * nBytes 

252 

253 # number of bits of the input 

254 inputType = im.dtype 

255 if inputType == numpy.uint8: 

256 nbi = 8 

257 elif inputType == numpy.uint16: 

258 nbi = 16 

259 elif inputType == numpy.uint32: 

260 nbi = 32 

261 elif inputType == numpy.uint64: 

262 nbi = 64 

263 else: 

264 raise ValueError(f"input image type should be unisgned integers ({inputType} given)") 

265 

266 return (float(2**nbo - 1) * (im) / float(2**nbi - 1)).astype("<u{}".format(nBytes)) 

267 

268 

269def singleShift(im, shift, axis, sub=0): 

270 """ 

271 This function shift the image and replace the border by an substitution value. 

272 

273 It uses ``numpy.roll``. 

274 

275 Parameters 

276 ----------- 

277 im : array 

278 The input to shift. 

279 shift : int 

280 The number of places by which elements are shifted (from numpy.roll). 

281 Default: 1 

282 axis : int 

283 The axis along which elements are shifted (from numpy.rool). 

284 sub : foat, default=0 

285 The substitution value of the border 

286 

287 Returns 

288 ------- 

289 array : 

290 The shifted image. 

291 

292 """ 

293 

294 # Step 1: Cyclic permutation on im 

295 im = numpy.roll(im, shift, axis=axis) 

296 

297 # Step 2: get image dimension 

298 dim = len(im.shape) 

299 

300 # Step 3: modify the boundary with replacement value 

301 if dim == 2: # if 2D image 

302 if shift == 1 and axis == 0: 

303 im[0, :] = sub 

304 elif shift == -1 and axis == 0: 

305 im[-1, :] = sub 

306 elif shift == 1 and axis == 1: 

307 im[:, 0] = sub 

308 elif shift == -1 and axis == 1: 

309 im[:, -1] = sub 

310 elif dim == 3: # if 3D image 

311 if shift >= 1 and axis == 0: 

312 im[0:shift, :, :] = sub 

313 elif shift <= -1 and axis == 0: 

314 im[shift:, :, :] = sub 

315 elif shift >= 1 and axis == 1: 

316 im[:, 0:shift, :] = sub 

317 elif shift <= -1 and axis == 1: 

318 im[:, shift:, :] = sub 

319 elif shift >= 1 and axis == 2: 

320 im[:, :, 0:shift] = sub 

321 elif shift <= -1 and axis == 2: 

322 im[:, :, shift:] = sub 

323 else: 

324 print("spam.helpers.imageManipulation.singleShift: dim={}. Should be 2 or 3.".format(dim)) 

325 print("exit function.") 

326 return -1 

327 

328 return im 

329 

330 

331def multipleShifts(im, shifts, sub=0): 

332 """ 

333 This function call ``singleShift`` multiple times. 

334 

335 Parameters 

336 ---------- 

337 im : array 

338 The input to shift. 

339 shifts : [int, int, int] 

340 Defines the number of shifts to apply in every axis. 

341 

342 .. code-block:: text 

343 

344 shift = [s_x, s_y, s_z] applies a shift of: 

345 . s_x on axis 0 

346 . s_y on axis 1 

347 . s_z on axis 2 

348 

349 sub : float, default=0 

350 The substitution value of the border 

351 

352 Returns 

353 ------- 

354 array : 

355 The shifted image. 

356 

357 """ 

358 

359 # loop over the n axis 

360 for i in range(len(shifts)): 

361 # if the value of the shift is not 0 on axis i 

362 if shifts[i]: 

363 # we call singleShift (only once) 

364 im = singleShift(im, shift=shifts[i], axis=i) 

365 

366 return im 

367 

368 

369def _binarisation(im, threshold=0.0, boolean=False, op=">", mask=None): 

370 """ 

371 This function binarise an input image according to a given threshold 

372 

373 It has an option to apply a mask to the binarized image to ignore the 

374 outside of the sample/specimen 

375 

376 Parameters 

377 ----------- 

378 im: array 

379 The image to binarise. 

380 

381 threshold : float 

382 the input limit value for binarization 

383 

384 boolean : bool 

385 Changes the output format and phase distribution (see output) 

386 

387 op : string, default='>' 

388 defines the thresholding operation 

389 

390 mask : array, default=None 

391 The mask of the input image: is 0 outside the boundary(specimen) and 1 inside 

392 

393 Returns 

394 -------- 

395 phases : array 

396 The repartition of phases resulting the binarisation. 

397 

398 For operator '>' it gives, if ``boolean=True``: 

399 

400 .. code-block:: text 

401 

402 0 - masked parts (where mask equals 0) 

403 0 - below threshold 

404 1 - above threshold 

405 

406 and if ``boolean=False`` 

407 

408 .. code-block:: text 

409 

410 0 - masked parts (where mask equals 0) 

411 1 - below threshold 

412 2 - above threshold 

413 """ 

414 

415 import operator 

416 

417 # Step 1: Get operator 

418 operation = { 

419 ">": operator.gt, 

420 "<": operator.lt, 

421 ">=": operator.ge, 

422 "<=": operator.le, 

423 "=": operator.eq, 

424 }.get(op) 

425 

426 # Step 2: binarisation 

427 phases = operation(im, threshold).astype("<u1") 

428 

429 # Step 3: rescaling if bool 

430 if not boolean: 

431 phases += 1 

432 

433 # Step 4: apply mask 

434 if mask is not None: 

435 phases = phases * mask 

436 

437 return phases 

438 

439 

440def slicePadded(im, startStop, createMask=False, padValue=0, verbose=False): 

441 """ 

442 Extract slice from im, padded with zeros, which is always of the dimensions asked (given from startStop) 

443 

444 Parameters 

445 ---------- 

446 im : 3D numpy array 

447 The image to be sliced 

448 

449 startStop : 6 component list of ints 

450 This array contains: 

451 [Zmin, Zmax, Ymin, Ymax, Xmin, Xmax] 

452 

453 createMask : bool, optional 

454 If True, return a padded slice, which is False when the slice falls outside im 

455 Default = False 

456 

457 Returns 

458 ------- 

459 imSliced : 3D numpy array 

460 The sliced image 

461 

462 mask : 3D numpy array of bools 

463 The 3D mask, only returned if createMask is True 

464 """ 

465 startStop = numpy.array(startStop).astype(int) 

466 

467 assert startStop[1] > startStop[0], "spam.helpers.slicePadded(): Zmax should be bigger than Zmin" 

468 assert startStop[3] > startStop[2], "spam.helpers.slicePadded(): Ymax should be bigger than Ymin" 

469 assert startStop[5] > startStop[4], "spam.helpers.slicePadded(): Xmax should be bigger than Xmin" 

470 

471 imSliced = ( 

472 numpy.zeros( 

473 ( 

474 startStop[1] - startStop[0], 

475 startStop[3] - startStop[2], 

476 startStop[5] - startStop[4], 

477 ), 

478 dtype=im.dtype, 

479 ) 

480 + padValue 

481 ) 

482 

483 start = numpy.array([startStop[0], startStop[2], startStop[4]]) 

484 stop = numpy.array([startStop[1], startStop[3], startStop[5]]) 

485 startOffset = numpy.array([max(0, start[0]), max(0, start[1]), max(0, start[2])]) 

486 stopOffset = numpy.array( 

487 [ 

488 min(im.shape[0], stop[0]), 

489 min(im.shape[1], stop[1]), 

490 min(im.shape[2], stop[2]), 

491 ] 

492 ) 

493 startLocal = startOffset - start 

494 stopLocal = startLocal + stopOffset - startOffset 

495 

496 # Check condition that we're asking for a slice of data wholly outside im 

497 # This means either that the stop values are all smaller than 0 

498 # OR the start are all bigger than im.shape 

499 if numpy.any(stop < numpy.array([0, 0, 0])) or numpy.any(start >= numpy.array(im.shape)): 

500 if verbose: 

501 print("spam.helpers.slicePadded(): The extracted padded slice doesn't not touch the image!") 

502 imSliced = imSliced.astype("<f4") 

503 imSliced *= numpy.nan 

504 if createMask: 

505 return imSliced, numpy.zeros_like(imSliced, dtype=bool) 

506 

507 else: 

508 imSliced[startLocal[0] : stopLocal[0], startLocal[1] : stopLocal[1], startLocal[2] : stopLocal[2],] = im[ 

509 startOffset[0] : stopOffset[0], 

510 startOffset[1] : stopOffset[1], 

511 startOffset[2] : stopOffset[2], 

512 ] 

513 if createMask: 

514 mask = numpy.zeros_like(imSliced, dtype=bool) 

515 mask[ 

516 startLocal[0] : stopLocal[0], 

517 startLocal[1] : stopLocal[1], 

518 startLocal[2] : stopLocal[2], 

519 ] = 1 

520 return imSliced, mask 

521 

522 return imSliced 

523 

524 

525def splitImage(im, divisions, margin, verbose=True): 

526 """ 

527 Divides the image in zDiv x yDiv x xDiv blocks, each block is padded 

528 with a margin. 

529 

530 Parameters 

531 ---------- 

532 im : 3D numpy array 

533 The image to be splitted 

534 

535 divisions : 3-component list of ints 

536 Desired number of blocks along Z, Y, X axes 

537 

538 margin : int 

539 Overlapping margin between each block. 

540 For applying a filter on subvolumes, it is recommended to use a margin of 1.5 times the filter diameter. 

541 For labelled data it is recommended that the margin is at least 1.5 times bigger than the particles largest axis 

542 

543 verbose : bool 

544 Print the parameters of the operations (number of blocks and margin) 

545 Default = True 

546 

547 Returns 

548 ------- 

549 Dictionary 

550 Dictionary with keys labelled acoording to the position of the block along each axis (e.g., 000, 001, 002,...) 

551 Each element (e.g., 001) within the dictionary carries the block origin and the resulting block, in that order 

552 

553 Note 

554 ---- 

555 This function should be used along `spam.helpers.imageManipulation.rebuildImage()` 

556 

557 """ 

558 zDiv, yDiv, xDiv = divisions 

559 

560 # Check if the slices can be made 

561 if zDiv >= im.shape[0]: 

562 print("spam.helpers.imageManipulation.splitImage: Incorrect number of slices for axis z") 

563 print("exit function.") 

564 return -1 

565 if yDiv >= im.shape[1]: 

566 print("spam.helpers.imageManipulation.splitImage: Incorrect number of slices for axis y") 

567 print("exit function.") 

568 return -1 

569 if xDiv >= im.shape[2]: 

570 print("spam.helpers.imageManipulation.splitImage: Incorrect number of slices for axis x") 

571 print("exit function.") 

572 return -1 

573 

574 # Check that margin is not greater than the slice 

575 if margin >= im.shape[0] / zDiv: 

576 print("spam.helpers.imageManipulation.splitImage: Margin is too big for z axis") 

577 print("exit function.") 

578 return -1 

579 if margin >= im.shape[1] / yDiv: 

580 print("spam.helpers.imageManipulation.splitImage: Margin is too big for y axis") 

581 print("exit function.") 

582 return -1 

583 if margin >= im.shape[2] / xDiv: 

584 print("spam.helpers.imageManipulation.splitImage: Margin is too big for x axis") 

585 print("exit function.") 

586 return -1 

587 # Print parameters if needed 

588 if verbose: 

589 print( 

590 "spam.helpers.imageManipulation.splitImage: Working with margin of ", 

591 margin, 

592 ". The total number of blocks is ", 

593 zDiv * yDiv * xDiv, 

594 ) 

595 

596 # Pad initial image with zeros on the edge 

597 imPad = numpy.pad(im, margin, mode="edge") 

598 

599 # Compute size of blocks 

600 zSize = int(im.shape[0] / zDiv) 

601 ySize = int(im.shape[1] / yDiv) 

602 xSize = int(im.shape[2] / xDiv) 

603 

604 # Create return dictionary 

605 output = {} 

606 

607 # Iterate through each block 

608 for zBlock in range(zDiv): 

609 for yBlock in range(yDiv): 

610 for xBlock in range(xDiv): 

611 # Get the origin of each block 

612 blockOrigin = numpy.array([zBlock * zSize, yBlock * ySize, xBlock * xSize]) 

613 blockSize = [zSize, ySize, xSize] 

614 # Check if the size needs to be changed to fit the image 

615 if zBlock == zDiv - 1 and im.shape[0] % zDiv != 0: 

616 blockSize[0] = zSize + (im.shape[0] % zDiv) 

617 if yBlock == yDiv - 1 and im.shape[1] % yDiv != 0: 

618 blockSize[1] = ySize + (im.shape[1] % yDiv) 

619 if xBlock == xDiv - 1 and im.shape[2] % xDiv != 0: 

620 blockSize[2] = xSize + (im.shape[2] % xDiv) 

621 

622 # Generate block with the margin on all sides 

623 imBlock = crop( 

624 imPad, 

625 ( 

626 blockSize[0] + 2 * margin, 

627 blockSize[1] + 2 * margin, 

628 blockSize[2] + 2 * margin, 

629 ), 

630 boxOrigin=(blockOrigin[0], blockOrigin[1], blockOrigin[2]), 

631 ) 

632 # Save the results 

633 output.update({str(zBlock) + str(yBlock) + str(xBlock): [blockOrigin, imBlock]}) 

634 # Save the margin 

635 output.update({"margin": margin}) 

636 

637 # Return 

638 return output 

639 

640 

641def rebuildImage(listBlocks, listCoordinates, margin, mode, keepLabels=False, verbose=True): 

642 """ 

643 Rebuilds splitted image from `spam.helpers.imageManipulation.splitImage()`. 

644 

645 Parameters 

646 ---------- 

647 listBlocks : list 

648 List of the 3D blocks that will form the re-built the image. 

649 Note: The order of listBlocks should be equivalent to the order of listCoordinates 

650 

651 listCoordinates : list 

652 List of the origin coordinates of each block. (Usually taken from `spam.helpers.imageManipulation.splitImage()`) 

653 Note: The order of listCoordinates should be equivalent to the order of listBlocks 

654 

655 margin : integer 

656 Value of the margin used for the images. (Usually taken from `spam.helpers.imageManipulation.splitImage()`) 

657 

658 mode : string 

659 'grey' : re-builds 3D greyscale arrays 

660 'label' : re-builds 3D labelled arrays 

661 

662 keepLabels : bool 

663 Do we need to want to keep the current labels from the blocks, or create a new one? 

664 Default = False 

665 

666 verbose : bool 

667 Print the evolution of the operation 

668 Default = True 

669 

670 Returns 

671 ------- 

672 imBuild : 3D numpy array 

673 Re-built image without the margins 

674 

675 Note 

676 ---- 

677 This function should be used along with `spam.helpers.imageManipulation.splitImage()` 

678 

679 """ 

680 

681 # Checking if listBlocks and listCoordinates have the same length 

682 if len(listBlocks) != len(listCoordinates): 

683 print("spam.helpers.imageManipulation.splitImage: listBlocks and listCoordinates must have the same length") 

684 return -1 

685 

686 # Transform listCoordinates into array 

687 arrayCoord = numpy.asarray(listCoordinates) 

688 

689 # Checking if all the origin coordinates are different 

690 _, counts = numpy.unique(arrayCoord, axis=0, return_counts=True) 

691 if len(counts) != len(arrayCoord): 

692 print("spam.helpers.imageManipulation.splitImage: coordinates in listCoordinates must be all different") 

693 return -1 

694 

695 if verbose: 

696 # Create progressbar 

697 widgets = [ 

698 progressbar.FormatLabel(""), 

699 " ", 

700 progressbar.Bar(), 

701 " ", 

702 progressbar.AdaptiveETA(), 

703 ] 

704 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(listCoordinates)) 

705 pbar.start() 

706 finishedBlocks = 0 

707 

708 if mode == "grey": 

709 

710 # Shape of the block opposite to the origine 

711 shapeLast = listBlocks[numpy.argmax(numpy.sum(arrayCoord, axis=1))].shape 

712 

713 # Tentative size of the final image 

714 zSize = numpy.amax(arrayCoord[:, 0]) + shapeLast[0] - 2 * margin 

715 ySize = numpy.amax(arrayCoord[:, 1]) + shapeLast[1] - 2 * margin 

716 xSize = numpy.amax(arrayCoord[:, 2]) + shapeLast[2] - 2 * margin 

717 

718 # Initialising rebuild image 

719 imBuild = numpy.zeros((zSize, ySize, xSize)) 

720 

721 # Loop on the length to lists, so to replace zeros in imBuild with the actual values at the right position 

722 for i in range(len(listCoordinates)): 

723 origin = listCoordinates[i] 

724 # GP 01/03/2022: Changing to ease the charge on the memory 

725 blockPad = listBlocks.pop(0) 

726 # blockPad = listBlocks[i] 

727 

728 if margin == 0: 

729 imBuild[ 

730 origin[0] : origin[0] + blockPad.shape[0] - 2 * margin, 

731 origin[1] : origin[1] + blockPad.shape[1] - 2 * margin, 

732 origin[2] : origin[2] + blockPad.shape[2] - 2 * margin, 

733 ] = blockPad 

734 else: 

735 imBuild[ 

736 origin[0] : origin[0] + blockPad.shape[0] - 2 * margin, 

737 origin[1] : origin[1] + blockPad.shape[1] - 2 * margin, 

738 origin[2] : origin[2] + blockPad.shape[2] - 2 * margin, 

739 ] = blockPad[margin:-margin, margin:-margin, margin:-margin] 

740 if verbose: 

741 # Update the progressbar 

742 finishedBlocks += 1 

743 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedBlocks, len(listCoordinates))) 

744 pbar.update(finishedBlocks) 

745 print("\n\t") 

746 return imBuild 

747 

748 if mode == "label": 

749 

750 # Shape of the block opposite to the origine 

751 shapeLast = listBlocks[numpy.argmax(numpy.sum(arrayCoord, axis=1))].shape 

752 

753 # Size of the final image 

754 zSize = numpy.amax(arrayCoord[:, 0]) + shapeLast[0] - 2 * margin 

755 ySize = numpy.amax(arrayCoord[:, 1]) + shapeLast[1] - 2 * margin 

756 xSize = numpy.amax(arrayCoord[:, 2]) + shapeLast[2] - 2 * margin 

757 

758 # Initialising rebuild image 

759 imBuild = numpy.zeros((zSize, ySize, xSize)) 

760 # Loop over the lists 

761 for i in range(len(listCoordinates)): 

762 # Get the origin of the block 

763 origin = listCoordinates[i] 

764 # Get the block 

765 # GP 01/03/2022: Changing to ease the charge on the memory 

766 # block = listBlocks[i] 

767 block = listBlocks.pop(0) 

768 # Compute the bounding boxes 

769 boundingBoxes = spam.label.boundingBoxes(block) 

770 # Compute centre of mass 

771 centresOfMass = spam.label.centresOfMass(block, boundingBoxes=boundingBoxes) 

772 

773 # List for classifying the labels 

774 inside = [] 

775 outside = [] 

776 partial = [] 

777 

778 # Check if each label is inside the true block - i.e., it is inside the block without the margin 

779 for j in range(1, len(boundingBoxes), 1): 

780 # Get the box 

781 box = boundingBoxes[j] 

782 # Check if the origin is inside the true block 

783 checkOrigin = margin < box[0] < block.shape[0] - margin and margin < box[2] < block.shape[1] - margin and margin < box[4] < block.shape[2] - margin 

784 # Check if the coordinate opposite to the origin is inside the true block 

785 checkOpp = margin < box[1] < block.shape[0] - margin and margin < box[3] < block.shape[1] - margin and margin < box[5] < block.shape[2] - margin 

786 if checkOrigin and checkOpp: 

787 # Both points are inside 

788 inside.append(j) 

789 elif checkOrigin or checkOpp: 

790 # Only one is inside 

791 partial.append(j) 

792 else: 

793 # Both are outside 

794 outside.append(j) 

795 

796 # Create true block array, keep only particles fully inside 

797 trueBlock = spam.label.removeLabels(block, partial + outside) 

798 

799 if not keepLabels: 

800 # Check that we have labels inside the block 

801 if len(numpy.unique(trueBlock)) > 1: 

802 # Make labels consecutive 

803 trueBlock = spam.label.makeLabelsSequential(trueBlock) 

804 trueBlock = numpy.where(trueBlock != 0, trueBlock + numpy.max(imBuild), trueBlock) 

805 

806 # IV (04-03-21): Info needed to avoid chopping grains that would be considered as outside particles 

807 imBuildSubSet = imBuild[ 

808 origin[0] : origin[0] + trueBlock.shape[0] - 2 * margin, 

809 origin[1] : origin[1] + trueBlock.shape[1] - 2 * margin, 

810 origin[2] : origin[2] + trueBlock.shape[2] - 2 * margin, 

811 ] 

812 

813 # Add the true block preserving previous info 

814 imBuild[origin[0] : origin[0] + trueBlock.shape[0] - 2 * margin, origin[1] : origin[1] + trueBlock.shape[1] - 2 * margin, origin[2] : origin[2] + trueBlock.shape[2] - 2 * margin,] = ( 

815 trueBlock[margin:-margin, margin:-margin, margin:-margin] + imBuildSubSet 

816 ) 

817 

818 # Get current maximum label 

819 tempMax = numpy.max(imBuild) 

820 # Label counter 

821 labCounter = 1 

822 # Solve the labels inside partial list 

823 if len(partial) > 0: 

824 # Iterate trough each label 

825 for label in partial: 

826 # Get the bounding box 

827 box = boundingBoxes[label] 

828 # Get subset 

829 labelSubset = spam.label.getLabel( 

830 block, 

831 label, 

832 boundingBoxes=boundingBoxes, 

833 centresOfMass=centresOfMass, 

834 ) 

835 labelSubvol = labelSubset["subvol"].astype(int) 

836 labelSubvol = numpy.where(labelSubvol != 0, labCounter + tempMax, labelSubvol) 

837 # IV+GP implementation - Just put the box there, overwrite whatever you find there 

838 # Get the subset from imBuild 

839 imBuildSubset = imBuild[ 

840 origin[0] + box[0] - margin : origin[0] + box[1] + 1 - margin, 

841 origin[1] + box[2] - margin : origin[1] + box[3] + 1 - margin, 

842 origin[2] + box[4] - margin : origin[2] + box[5] + 1 - margin, 

843 ] 

844 # Change 

845 imBuildSubset = numpy.where(labelSubvol != 0, labelSubvol, imBuildSubset) 

846 # Put back 

847 imBuild[ 

848 origin[0] + box[0] - margin : origin[0] + box[1] + 1 - margin, 

849 origin[1] + box[2] - margin : origin[1] + box[3] + 1 - margin, 

850 origin[2] + box[4] - margin : origin[2] + box[5] + 1 - margin, 

851 ] = imBuildSubset 

852 # Update label counter 

853 labCounter += 1 

854 

855 if verbose: 

856 # Update the progressbar 

857 finishedBlocks += 1 

858 widgets[0] = progressbar.FormatLabel("{}/{} ".format(finishedBlocks, len(listCoordinates))) 

859 pbar.update(finishedBlocks)