Coverage for openhcs/processing/backends/processors/cupy_processor.py: 15.5%

1""" 

2CuPy Image Processor Implementation 

3 

4This module implements the ImageProcessorInterface using CuPy as the backend. 

5It leverages GPU acceleration for image processing operations. 

6 

7Doctrinal Clauses: 

8- Clause 3 — Declarative Primacy: All functions are pure and stateless 

9- Clause 88 — No Inferred Capabilities: Explicit CuPy dependency 

10- Clause 106-A — Declared Memory Types: All methods specify CuPy arrays 

11""" 

12from __future__ import annotations 

13 

14import logging 

15import os 

16from typing import Any, List, Optional, Tuple 

17 

18from openhcs.core.memory.decorators import cupy as cupy_func 

19from openhcs.core.utils import optional_import 

20 

21logger = logging.getLogger(__name__) 

22 

23# Check if we're in subprocess runner mode and should skip GPU imports 

24if os.getenv('OPENHCS_SUBPROCESS_NO_GPU') == '1': 24 ↛ 26line 24 didn't jump to line 26 because the condition on line 24 was never true

25 # Subprocess runner mode - skip GPU imports 

26 cp = None 

27 ndimage = None 

28 cucim_filters = None 

29 logger.info("Subprocess runner mode - skipping cupy import") 

30else: 

31 # Normal mode - import CuPy as an optional dependency 

32 cp = optional_import("cupy") 

33 ndimage = None 

34 if cp is not None: 34 ↛ 40line 34 didn't jump to line 40 because the condition on line 34 was always true

35 cupyx_scipy = optional_import("cupyx.scipy") 

36 if cupyx_scipy is not None: 36 ↛ 40line 36 didn't jump to line 40 because the condition on line 36 was always true

37 ndimage = cupyx_scipy.ndimage 

38 

39 # Import CuCIM for edge detection 

40 cucim_filters = optional_import("cucim.skimage.filters") 

41 

42logger = logging.getLogger(__name__) 

43 

44 

45@cupy_func 

46def create_linear_weight_mask(height: int, width: int, margin_ratio: float = 0.1) -> "cp.ndarray": 

47 """ 

48 Create a 2D weight mask that linearly ramps from 0 at the edges to 1 in the center. 

49 

50 Args: 

51 height: Height of the mask 

52 width: Width of the mask 

53 margin_ratio: Ratio of the margin to the image size 

54 

55 Returns: 

56 2D CuPy weight mask of shape (height, width) 

57 """ 

58 # The compiler will ensure this function is only called when CuPy is available 

59 # No need to check for CuPy availability here 

60 

61 margin_y = int(cp.floor(height * margin_ratio)) 

62 margin_x = int(cp.floor(width * margin_ratio)) 

63 

64 weight_y = cp.ones(height, dtype=cp.float32) 

65 if margin_y > 0: 

66 ramp_top = cp.linspace(0, 1, margin_y, endpoint=False) 

67 ramp_bottom = cp.linspace(1, 0, margin_y, endpoint=False) 

68 weight_y[:margin_y] = ramp_top 

69 weight_y[-margin_y:] = ramp_bottom 

70 

71 weight_x = cp.ones(width, dtype=cp.float32) 

72 if margin_x > 0: 

73 ramp_left = cp.linspace(0, 1, margin_x, endpoint=False) 

74 ramp_right = cp.linspace(1, 0, margin_x, endpoint=False) 

75 weight_x[:margin_x] = ramp_left 

76 weight_x[-margin_x:] = ramp_right 

77 

78 # Create 2D weight mask 

79 weight_mask = cp.outer(weight_y, weight_x) 

80 

81 return weight_mask 

82 

83 

84def _validate_3d_array(array: Any, name: str = "input") -> None: 

85 """ 

86 Validate that the input is a 3D CuPy array. 

87 

88 Args: 

89 array: Array to validate 

90 name: Name of the array for error messages 

91 

92 Raises: 

93 TypeError: If the array is not a CuPy array 

94 ValueError: If the array is not 3D 

95 ImportError: If CuPy is not available 

96 """ 

97 # The compiler will ensure this function is only called when CuPy is available 

98 # No need to check for CuPy availability here 

99 

100 if not isinstance(array, cp.ndarray): 

101 raise TypeError(f"{name} must be a CuPy array, got {type(array)}. " 

102 f"No automatic conversion is performed to maintain explicit contracts.") 

103 

104 if array.ndim != 3: 

105 raise ValueError(f"{name} must be a 3D array, got {array.ndim}D") 

106 

107@cupy_func 

108def sharpen(image: "cp.ndarray", radius: float = 1.0, amount: float = 1.0) -> "cp.ndarray": 

109 """ 

110 Sharpen a 3D image using unsharp masking - GPU PARALLELIZED. 

111 

112 This applies sharpening to each Z-slice independently using vectorized operations 

113 for maximum GPU utilization. Normalization and rescaling are fully parallelized. 

114 

115 Args: 

116 image: 3D CuPy array of shape (Z, Y, X) 

117 radius: Radius of Gaussian blur 

118 amount: Sharpening strength 

119 

120 Returns: 

121 Sharpened 3D CuPy array of shape (Z, Y, X) 

122 """ 

123 _validate_3d_array(image) 

124 

125 # Store original dtype 

126 dtype = image.dtype 

127 

128 # Convert to float32 for processing 

129 image_float = image.astype(cp.float32) 

130 

131 # Vectorized per-slice normalization - GPU PARALLELIZED 

132 # Get max value per slice: shape (Z, 1, 1) for broadcasting 

133 max_per_slice = cp.max(image_float, axis=(1, 2), keepdims=True) 

134 image_norm = image_float / max_per_slice 

135 

136 # Apply 3D Gaussian blur with sigma_z=0 for slice-wise processing - GPU PARALLELIZED 

137 # This processes all slices simultaneously while keeping Z-slices independent 

138 blurred = ndimage.gaussian_filter(image_norm, sigma=(0, radius, radius)) 

139 

140 # Apply unsharp mask: original + amount * (original - blurred) - GPU PARALLELIZED 

141 sharpened = image_norm + amount * (image_norm - blurred) 

142 

143 # Clip to valid range - GPU PARALLELIZED 

144 sharpened = cp.clip(sharpened, 0, 1.0) 

145 

146 # Vectorized rescaling back to original range - GPU PARALLELIZED 

147 min_per_slice = cp.min(sharpened, axis=(1, 2), keepdims=True) 

148 max_per_slice = cp.max(sharpened, axis=(1, 2), keepdims=True) 

149 

150 # Avoid division by zero using broadcasting 

151 range_per_slice = max_per_slice - min_per_slice 

152 # Only rescale slices where max > min 

153 valid_range = range_per_slice > 0 

154 sharpened_rescaled = cp.where( 

155 valid_range, 

156 (sharpened - min_per_slice) * 65535 / range_per_slice, 

157 sharpened * 65535 

158 ) 

159 

160 # Convert back to original dtype 

161 return sharpened_rescaled.astype(dtype) 

162 

163@cupy_func 

164def percentile_normalize( 

165 image: "cp.ndarray", 

166 low_percentile: float = 1.0, 

167 high_percentile: float = 99.0, 

168 target_min: float = None, 

169 target_max: float = None, 

170 preserve_dtype: bool = True 

171) -> "cp.ndarray": 

172 """ 

173 Normalize a 3D image using percentile-based contrast stretching. 

174 

175 This applies normalization to each Z-slice independently using slice-by-slice 

176 processing for algorithmic reasons (each slice needs different percentile values). 

177 Each slice operation is GPU parallelized across Y,X dimensions. 

178 

179 Args: 

180 image: 3D CuPy array of shape (Z, Y, X) 

181 low_percentile: Lower percentile (0-100) 

182 high_percentile: Upper percentile (0-100) 

183 target_min: Target minimum value (auto-detected from dtype if None) 

184 target_max: Target maximum value (auto-detected from dtype if None) 

185 preserve_dtype: If True, output same dtype as input. If False, use target range. 

186 

187 Returns: 

188 Normalized 3D CuPy array of shape (Z, Y, X) with same or specified dtype 

189 """ 

190 _validate_3d_array(image) 

191 

192 # Import shared utilities 

193 from .percentile_utils import resolve_target_range, slice_percentile_normalize_core 

194 

195 # Auto-detect target range based on input dtype if not specified 

196 target_min, target_max = resolve_target_range(image.dtype, target_min, target_max) 

197 

198 # Use shared core logic with CuPy-specific functions 

199 return slice_percentile_normalize_core( 

200 image=image, 

201 low_percentile=low_percentile, 

202 high_percentile=high_percentile, 

203 target_min=target_min, 

204 target_max=target_max, 

205 percentile_func=cp.percentile, 

206 clip_func=cp.clip, 

207 ones_like_func=cp.ones_like, 

208 zeros_like_func=lambda arr, dtype=None: cp.zeros_like(arr, dtype=dtype or cp.float32), 

209 preserve_dtype=preserve_dtype 

210 ) 

211 

212@cupy_func 

213def stack_percentile_normalize( 

214 stack: "cp.ndarray", 

215 low_percentile: float = 1.0, 

216 high_percentile: float = 99.0, 

217 target_min: float = None, 

218 target_max: float = None, 

219 preserve_dtype: bool = True 

220) -> "cp.ndarray": 

221 """ 

222 Normalize a stack using global percentile-based contrast stretching. 

223 

224 This ensures consistent normalization across all Z-slices by computing 

225 global percentiles across the entire stack. 

226 

227 Args: 

228 stack: 3D CuPy array of shape (Z, Y, X) 

229 low_percentile: Lower percentile (0-100) 

230 high_percentile: Upper percentile (0-100) 

231 target_min: Target minimum value (auto-detected from dtype if None) 

232 target_max: Target maximum value (auto-detected from dtype if None) 

233 preserve_dtype: If True, output same dtype as input. If False, use target range. 

234 

235 Returns: 

236 Normalized 3D CuPy array of shape (Z, Y, X) with same or specified dtype 

237 """ 

238 _validate_3d_array(stack) 

239 

240 # Import shared utilities 

241 from .percentile_utils import resolve_target_range, percentile_normalize_core 

242 

243 # Auto-detect target range based on input dtype if not specified 

244 target_min, target_max = resolve_target_range(stack.dtype, target_min, target_max) 

245 

246 # Use shared core logic with CuPy-specific functions 

247 return percentile_normalize_core( 

248 stack=stack, 

249 low_percentile=low_percentile, 

250 high_percentile=high_percentile, 

251 target_min=target_min, 

252 target_max=target_max, 

253 percentile_func=lambda arr, pct: cp.percentile(arr, pct), 

254 clip_func=cp.clip, 

255 ones_like_func=cp.ones_like, 

256 preserve_dtype=preserve_dtype 

257 ) 

258 

259@cupy_func 

260def create_composite( 

261 stack: "cp.ndarray", weights: Optional[List[float]] = None 

262) -> "cp.ndarray": 

263 """ 

264 Create a composite image from a 3D stack where each slice is a channel. 

265 

266 Args: 

267 stack: 3D CuPy array of shape (N, Y, X) where N is number of channel slices 

268 weights: List of weights for each slice. If None, equal weights are used. 

269 

270 Returns: 

271 Composite 3D CuPy array of shape (1, Y, X) 

272 """ 

273 # 🔄 MEMORY CONVERSION LOGGING: Log what we actually received 

274 import logging 

275 logger = logging.getLogger(__name__) 

276 logger.info(f"🔄 CREATE_COMPOSITE: Called with stack type: {type(stack)}, shape: {getattr(stack, 'shape', 'no shape')}") 

277 logger.info(f"🔄 CREATE_COMPOSITE: weights parameter - type: {type(weights)}, value: {weights}") 

278 

279 # Validate input is 3D array 

280 if not hasattr(stack, 'shape') or len(stack.shape) != 3: 

281 raise TypeError(f"stack must be a 3D CuPy array, got shape: {getattr(stack, 'shape', 'no shape')}") 

282 

283 n_slices, height, width = stack.shape 

284 

285 # Default weights if none provided 

286 if weights is None: 

287 # Equal weights for all slices 

288 weights = [1.0 / n_slices] * n_slices 

289 logger.info(f"🔄 CREATE_COMPOSITE: Using default equal weights: {weights}") 

290 elif isinstance(weights, (list, tuple)): 

291 # Convert tuple to list if needed 

292 weights = list(weights) 

293 logger.info(f"🔄 CREATE_COMPOSITE: Using provided weights: {weights}") 

294 if len(weights) != n_slices: 

295 raise ValueError(f"Number of weights ({len(weights)}) must match number of slices ({n_slices})") 

296 else: 

297 # Log the problematic type and value for debugging 

298 logger.error(f"🔄 CREATE_COMPOSITE: Invalid weights type - expected list/tuple/None, got {type(weights)}: {weights}") 

299 raise TypeError(f"weights must be a list of values or None, got {type(weights)}: {weights}") 

300 

301 # Normalize weights to sum to 1 

302 weight_sum = sum(weights) 

303 if weight_sum == 0: 

304 raise ValueError("Sum of weights cannot be zero") 

305 normalized_weights = [w / weight_sum for w in weights] 

306 

307 # Convert weights to CuPy array for efficient computation 

308 # CRITICAL: Use float32 for weights to preserve fractional values, not stack.dtype 

309 weights_array = cp.array(normalized_weights, dtype=cp.float32) 

310 

311 # Reshape weights for broadcasting: (N, 1, 1) to multiply with (N, Y, X) 

312 weights_array = weights_array.reshape(n_slices, 1, 1) 

313 

314 # Create composite by weighted sum along the first axis 

315 # Convert stack to float32 for computation to avoid precision loss 

316 stack_float = stack.astype(cp.float32) 

317 weighted_stack = stack_float * weights_array 

318 composite_slice = cp.sum(weighted_stack, axis=0, keepdims=True) # Keep as (1, Y, X) 

319 

320 # Convert back to original dtype 

321 composite_slice = composite_slice.astype(stack.dtype) 

322 

323 return composite_slice 

324 

325@cupy_func 

326def apply_mask(image: "cp.ndarray", mask: "cp.ndarray") -> "cp.ndarray": 

327 """ 

328 Apply a mask to a 3D image. 

329 

330 This applies the mask to each Z-slice independently if mask is 2D, 

331 or applies the 3D mask directly if mask is 3D. 

332 

333 Args: 

334 image: 3D CuPy array of shape (Z, Y, X) 

335 mask: 3D CuPy array of shape (Z, Y, X) or 2D CuPy array of shape (Y, X) 

336 

337 Returns: 

338 Masked 3D CuPy array of shape (Z, Y, X) 

339 """ 

340 _validate_3d_array(image) 

341 

342 # Handle 2D mask (apply to each Z-slice) 

343 if isinstance(mask, cp.ndarray) and mask.ndim == 2: 

344 if mask.shape != image.shape[1:]: 

345 raise ValueError( 

346 f"2D mask shape {mask.shape} doesn't match image slice shape {image.shape[1:]}" 

347 ) 

348 

349 # Apply 2D mask to all Z-slices simultaneously using broadcasting - GPU PARALLELIZED 

350 # Broadcasting mask from (Y, X) to (1, Y, X) allows vectorized operation across all slices 

351 mask_3d = mask[None, :, :] # Shape: (1, Y, X) 

352 result = image.astype(cp.float32) * mask_3d.astype(cp.float32) 

353 

354 return result.astype(image.dtype) 

355 

356 # Handle 3D mask 

357 if isinstance(mask, cp.ndarray) and mask.ndim == 3: 

358 if mask.shape != image.shape: 

359 raise ValueError( 

360 f"3D mask shape {mask.shape} doesn't match image shape {image.shape}" 

361 ) 

362 

363 # Apply 3D mask directly - MATCH NUMPY EXACTLY 

364 masked = image.astype(cp.float32) * mask.astype(cp.float32) 

365 return masked.astype(image.dtype) 

366 

367 # If we get here, the mask is neither 2D nor 3D CuPy array 

368 raise TypeError(f"mask must be a 2D or 3D CuPy array, got {type(mask)}") 

369 

370@cupy_func 

371def create_weight_mask(shape: Tuple[int, int], margin_ratio: float = 0.1) -> "cp.ndarray": 

372 """ 

373 Create a weight mask for blending images. 

374 

375 Args: 

376 shape: Shape of the mask (height, width) 

377 margin_ratio: Ratio of image size to use as margin 

378 

379 Returns: 

380 2D CuPy weight mask of shape (Y, X) 

381 """ 

382 if not isinstance(shape, tuple) or len(shape) != 2: 

383 raise TypeError("shape must be a tuple of (height, width)") 

384 

385 height, width = shape 

386 return create_linear_weight_mask(height, width, margin_ratio) 

387 

388@cupy_func 

389def max_projection(stack: "cp.ndarray") -> "cp.ndarray": 

390 """ 

391 Create a maximum intensity projection from a Z-stack. 

392 

393 Args: 

394 stack: 3D CuPy array of shape (Z, Y, X) 

395 

396 Returns: 

397 3D CuPy array of shape (1, Y, X) 

398 """ 

399 _validate_3d_array(stack) 

400 

401 # Create max projection 

402 projection_2d = cp.max(stack, axis=0) 

403 return projection_2d.reshape(1, projection_2d.shape[0], projection_2d.shape[1]) 

404 

405@cupy_func 

406def mean_projection(stack: "cp.ndarray") -> "cp.ndarray": 

407 """ 

408 Create a mean intensity projection from a Z-stack. 

409 

410 Args: 

411 stack: 3D CuPy array of shape (Z, Y, X) 

412 

413 Returns: 

414 3D CuPy array of shape (1, Y, X) 

415 """ 

416 _validate_3d_array(stack) 

417 

418 # Create mean projection 

419 projection_2d = cp.mean(stack, axis=0).astype(stack.dtype) 

420 return projection_2d.reshape(1, projection_2d.shape[0], projection_2d.shape[1]) 

421 

422@cupy_func 

423def spatial_bin_2d( 

424 stack: "cp.ndarray", 

425 bin_size: int = 2, 

426 method: str = "mean" 

427) -> "cp.ndarray": 

428 """ 

429 Apply 2D spatial binning to each slice in the stack - GPU accelerated. 

430 

431 Reduces spatial resolution by combining neighboring pixels in 2D blocks. 

432 Each slice is processed independently using efficient GPU operations. 

433 

434 Args: 

435 stack: 3D CuPy array of shape (Z, Y, X) 

436 bin_size: Size of the square binning kernel (e.g., 2 = 2x2 binning) 

437 method: Binning method - "mean", "sum", "max", or "min" 

438 

439 Returns: 

440 Binned 3D CuPy array of shape (Z, Y//bin_size, X//bin_size) 

441 """ 

442 _validate_3d_array(stack) 

443 

444 if bin_size <= 0: 

445 raise ValueError("bin_size must be positive") 

446 if method not in ["mean", "sum", "max", "min"]: 

447 raise ValueError("method must be one of: mean, sum, max, min") 

448 

449 z_slices, height, width = stack.shape 

450 

451 # Calculate output dimensions 

452 new_height = height // bin_size 

453 new_width = width // bin_size 

454 

455 if new_height == 0 or new_width == 0: 

456 raise ValueError(f"bin_size {bin_size} is too large for image dimensions {height}x{width}") 

457 

458 # Crop to make dimensions divisible by bin_size 

459 crop_height = new_height * bin_size 

460 crop_width = new_width * bin_size 

461 cropped_stack = stack[:, :crop_height, :crop_width] 

462 

463 # Reshape for binning: (Z, new_height, bin_size, new_width, bin_size) 

464 reshaped = cropped_stack.reshape(z_slices, new_height, bin_size, new_width, bin_size) 

465 

466 # Apply binning operation using CuPy functions 

467 if method == "mean": 

468 result = cp.mean(reshaped, axis=(2, 4)) 

469 elif method == "sum": 

470 result = cp.sum(reshaped, axis=(2, 4)) 

471 elif method == "max": 

472 result = cp.max(reshaped, axis=(2, 4)) 

473 elif method == "min": 

474 result = cp.min(reshaped, axis=(2, 4)) 

475 

476 return result.astype(stack.dtype) 

477 

478@cupy_func 

479def spatial_bin_3d( 

480 stack: "cp.ndarray", 

481 bin_size: int = 2, 

482 method: str = "mean" 

483) -> "cp.ndarray": 

484 """ 

485 Apply 3D spatial binning to the entire stack - GPU accelerated. 

486 

487 Reduces spatial resolution by combining neighboring voxels in 3D blocks 

488 using efficient GPU operations. 

489 

490 Args: 

491 stack: 3D CuPy array of shape (Z, Y, X) 

492 bin_size: Size of the cubic binning kernel (e.g., 2 = 2x2x2 binning) 

493 method: Binning method - "mean", "sum", "max", or "min" 

494 

495 Returns: 

496 Binned 3D CuPy array of shape (Z//bin_size, Y//bin_size, X//bin_size) 

497 """ 

498 _validate_3d_array(stack) 

499 

500 if bin_size <= 0: 

501 raise ValueError("bin_size must be positive") 

502 if method not in ["mean", "sum", "max", "min"]: 

503 raise ValueError("method must be one of: mean, sum, max, min") 

504 

505 depth, height, width = stack.shape 

506 

507 # Calculate output dimensions 

508 new_depth = depth // bin_size 

509 new_height = height // bin_size 

510 new_width = width // bin_size 

511 

512 if new_depth == 0 or new_height == 0 or new_width == 0: 

513 raise ValueError(f"bin_size {bin_size} is too large for stack dimensions {depth}x{height}x{width}") 

514 

515 # Crop to make dimensions divisible by bin_size 

516 crop_depth = new_depth * bin_size 

517 crop_height = new_height * bin_size 

518 crop_width = new_width * bin_size 

519 cropped_stack = stack[:crop_depth, :crop_height, :crop_width] 

520 

521 # Reshape for 3D binning: (new_depth, bin_size, new_height, bin_size, new_width, bin_size) 

522 reshaped = cropped_stack.reshape(new_depth, bin_size, new_height, bin_size, new_width, bin_size) 

523 

524 # Apply binning operation across the three bin_size dimensions using CuPy functions 

525 if method == "mean": 

526 result = cp.mean(reshaped, axis=(1, 3, 5)) 

527 elif method == "sum": 

528 result = cp.sum(reshaped, axis=(1, 3, 5)) 

529 elif method == "max": 

530 result = cp.max(reshaped, axis=(1, 3, 5)) 

531 elif method == "min": 

532 result = cp.min(reshaped, axis=(1, 3, 5)) 

533 

534 return result.astype(stack.dtype) 

535 

536@cupy_func 

537def stack_equalize_histogram( 

538 stack: "cp.ndarray", 

539 bins: int = 65536, 

540 range_min: float = 0.0, 

541 range_max: float = 65535.0 

542) -> "cp.ndarray": 

543 """ 

544 Apply histogram equalization to an entire stack. 

545 

546 This ensures consistent contrast enhancement across all Z-slices by 

547 computing a global histogram across the entire stack. 

548 

549 Args: 

550 stack: 3D CuPy array of shape (Z, Y, X) 

551 bins: Number of bins for histogram computation 

552 range_min: Minimum value for histogram range 

553 range_max: Maximum value for histogram range 

554 

555 Returns: 

556 Equalized 3D CuPy array of shape (Z, Y, X) 

557 """ 

558 _validate_3d_array(stack) 

559 

560 # MATCH NUMPY EXACTLY - Flatten the entire stack to compute the global histogram 

561 flat_stack = stack.flatten() 

562 

563 # Calculate the histogram and cumulative distribution function (CDF) - MATCH NUMPY EXACTLY 

564 hist, bin_edges = cp.histogram(flat_stack, bins=bins, range=(range_min, range_max)) 

565 cdf = hist.cumsum() 

566 

567 # Normalize the CDF to the range [0, 65535] - MATCH NUMPY EXACTLY 

568 # Avoid division by zero 

569 if cdf[-1] > 0: 

570 cdf = 65535 * cdf / cdf[-1] 

571 

572 # Use linear interpolation to map input values to equalized values - MATCH NUMPY EXACTLY 

573 equalized_stack = cp.interp(stack.flatten(), bin_edges[:-1], cdf).reshape(stack.shape) 

574 

575 # Convert to uint16 - MATCH NUMPY EXACTLY 

576 return equalized_stack.astype(cp.uint16) 

577 

578@cupy_func 

579def create_projection( 

580 stack: "cp.ndarray", method: str = "max_projection" 

581) -> "cp.ndarray": 

582 """ 

583 Create a projection from a stack using the specified method. 

584 

585 Args: 

586 stack: 3D CuPy array of shape (Z, Y, X) 

587 method: Projection method (max_projection, mean_projection) 

588 

589 Returns: 

590 3D CuPy array of shape (1, Y, X) 

591 """ 

592 _validate_3d_array(stack) 

593 

594 if method == "max_projection": 

595 return max_projection(stack) 

596 

597 if method == "mean_projection": 

598 return mean_projection(stack) 

599 

600 # FAIL FAST: No fallback projection methods 

601 raise ValueError(f"Unknown projection method: {method}. Valid methods: max_projection, mean_projection") 

602 

603 

604@cupy_func 

605def crop( 

606 input_image: "cp.ndarray", 

607 start_x: int = 0, 

608 start_y: int = 0, 

609 start_z: int = 0, 

610 width: int = 1, 

611 height: int = 1, 

612 depth: int = 1 

613) -> "cp.ndarray": 

614 """ 

615 Crop a given substack out of a given image stack. 

616 

617 Equivalent to pyclesperanto.crop() but using CuPy operations. 

618 

619 Parameters 

620 ---------- 

621 input_image: cp.ndarray 

622 Input 3D image to process of shape (Z, Y, X) 

623 start_x: int (= 0) 

624 Starting index coordinate x 

625 start_y: int (= 0) 

626 Starting index coordinate y 

627 start_z: int (= 0) 

628 Starting index coordinate z 

629 width: int (= 1) 

630 Width size of the region to crop 

631 height: int (= 1) 

632 Height size of the region to crop 

633 depth: int (= 1) 

634 Depth size of the region to crop 

635 

636 Returns 

637 ------- 

638 cp.ndarray 

639 Cropped 3D array of shape (depth, height, width) 

640 """ 

641 _validate_3d_array(input_image) 

642 

643 # Validate crop parameters 

644 if width <= 0 or height <= 0 or depth <= 0: 

645 raise ValueError(f"Crop dimensions must be positive: width={width}, height={height}, depth={depth}") 

646 

647 if start_x < 0 or start_y < 0 or start_z < 0: 

648 raise ValueError(f"Start coordinates must be non-negative: start_x={start_x}, start_y={start_y}, start_z={start_z}") 

649 

650 # Get input dimensions 

651 input_depth, input_height, input_width = input_image.shape 

652 

653 # Calculate end coordinates 

654 end_x = start_x + width 

655 end_y = start_y + height 

656 end_z = start_z + depth 

657 

658 # Validate bounds 

659 if end_x > input_width or end_y > input_height or end_z > input_depth: 

660 raise ValueError( 

661 f"Crop region extends beyond image bounds. " 

662 f"Image shape: {input_image.shape}, " 

663 f"Crop region: ({start_z}:{end_z}, {start_y}:{end_y}, {start_x}:{end_x})" 

664 ) 

665 

666 # Perform the crop using CuPy slicing 

667 cropped = input_image[start_z:end_z, start_y:end_y, start_x:end_x] 

668 

669 return cropped 

670 

671def _create_disk_cupy(radius: int) -> "cp.ndarray": 

672 """Create a disk structuring element using CuPy - MATCH NUMPY EXACTLY""" 

673 y, x = cp.ogrid[-radius:radius+1, -radius:radius+1] 

674 mask = x*x + y*y <= radius*radius 

675 return mask.astype(cp.uint8) 

676 

677def _resize_cupy_better_match(image: "cp.ndarray", output_shape: tuple, anti_aliasing: bool = True, preserve_range: bool = True) -> "cp.ndarray": 

678 """ 

679 Resize image using CuPy to better match scikit-image's transform.resize behavior. 

680 

681 Key differences from original: 

682 1. Better anti-aliasing sigma calculation 

683 2. Proper preserve_range handling without rescaling 

684 3. More accurate zoom parameters 

685 """ 

686 from cupyx.scipy import ndimage as cupy_ndimage 

687 

688 # Calculate zoom factors 

689 zoom_factors = [output_shape[i] / image.shape[i] for i in range(len(output_shape))] 

690 

691 # Convert to float32 for processing 

692 image_float = image.astype(cp.float32) 

693 

694 # Apply anti-aliasing if downsampling 

695 if anti_aliasing and any(z < 1.0 for z in zoom_factors): 

696 # Use scikit-image's sigma calculation: sigma = (1/zoom - 1) / 2 

697 # But ensure minimum sigma and handle edge cases 

698 sigma = [] 

699 for z in zoom_factors: 

700 if z < 1.0: 

701 s = (1.0/z - 1.0) / 2.0 

702 sigma.append(max(s, 0.5)) # Minimum sigma for stability 

703 else: 

704 sigma.append(0.0) 

705 

706 # Apply Gaussian smoothing before downsampling 

707 if any(s > 0 for s in sigma): 

708 image_float = cupy_ndimage.gaussian_filter(image_float, sigma) 

709 

710 # Perform zoom with bilinear interpolation (order=1) 

711 resized = cupy_ndimage.zoom(image_float, zoom_factors, order=1) 

712 

713 # Handle preserve_range properly - don't rescale, just maintain dtype range 

714 if preserve_range: 

715 # Clip to valid range for the original dtype 

716 if image.dtype == cp.uint16: 

717 resized = cp.clip(resized, 0, 65535) 

718 elif image.dtype == cp.uint8: 

719 resized = cp.clip(resized, 0, 255) 

720 # For other dtypes, keep as-is 

721 

722 return resized 

723 

724@cupy_func 

725def tophat( 

726 image: "cp.ndarray", 

727 selem_radius: int = 50, 

728 downsample_factor: int = 4, 

729 downsample_anti_aliasing: bool = True, 

730 upsample_anti_aliasing: bool = False 

731) -> "cp.ndarray": 

732 """ 

733 Apply white top-hat filter to a 3D image for background removal. 

734 

735 This applies the filter to each Z-slice independently using slice-by-slice 

736 processing for algorithmic reasons (complex multi-step processing with 

737 slice-specific intermediate results). Each slice operation is GPU parallelized. 

738 

739 Args: 

740 image: 3D CuPy array of shape (Z, Y, X) 

741 selem_radius: Radius of the structuring element disk 

742 downsample_factor: Factor by which to downsample the image for processing 

743 downsample_anti_aliasing: Whether to use anti-aliasing when downsampling 

744 upsample_anti_aliasing: Whether to use anti-aliasing when upsampling 

745 

746 Returns: 

747 Filtered 3D CuPy array of shape (Z, Y, X) 

748 """ 

749 _validate_3d_array(image) 

750 

751 # Process each Z-slice independently - IMPROVED MATCH TO NUMPY 

752 # NOTE: For loop is ALGORITHMICALLY NECESSARY here due to complex multi-step 

753 # processing (downsample, morphology, upsample, subtract) with slice-specific 

754 # intermediate results. Vectorization would require significant algorithm restructuring. 

755 result = cp.zeros_like(image) 

756 

757 for z in range(image.shape[0]): 

758 # Store original data type - MATCH NUMPY EXACTLY 

759 input_dtype = image[z].dtype 

760 

761 # 1) Downsample - IMPROVED MATCH TO NUMPY 

762 target_shape = (image[z].shape[0]//downsample_factor, image[z].shape[1]//downsample_factor) 

763 image_small = _resize_cupy_better_match( 

764 image[z], 

765 target_shape, 

766 anti_aliasing=downsample_anti_aliasing, 

767 preserve_range=True 

768 ) 

769 

770 # 2) Build structuring element for the smaller image - MATCH NUMPY EXACTLY 

771 selem_small = _create_disk_cupy(selem_radius // downsample_factor) 

772 

773 # 3) White top-hat on the smaller image - MATCH NUMPY EXACTLY 

774 tophat_small = ndimage.white_tophat(image_small, structure=selem_small) 

775 

776 # 4) Upscale background to original size - IMPROVED MATCH TO NUMPY 

777 background_small = image_small - tophat_small 

778 background_large = _resize_cupy_better_match( 

779 background_small, 

780 image[z].shape, 

781 anti_aliasing=upsample_anti_aliasing, 

782 preserve_range=True 

783 ) 

784 

785 # 5) Subtract background and clip negative values - MATCH NUMPY EXACTLY 

786 slice_result = cp.maximum(image[z] - background_large, 0) 

787 

788 # 6) Convert back to original data type - MATCH NUMPY EXACTLY 

789 result[z] = slice_result.astype(input_dtype) 

790 

791 return result 

792 

793# Lazy initialization of ElementwiseKernel to avoid import errors 

794_sobel_2d_parallel = None 

795 

796def _get_sobel_2d_kernel(): 

797 """Get or create the Sobel 2D parallel kernel.""" 

798 global _sobel_2d_parallel 

799 if _sobel_2d_parallel is None: 

800 if cp is None: 

801 raise ImportError("CuPy is required for GPU Sobel operations") 

802 

803 _sobel_2d_parallel = cp.ElementwiseKernel( 

804 'raw T input, int32 Y, int32 X, int32 mode, T cval', 

805 'T output', 

806 ''' 

807 // Calculate which slice, row, col this thread handles 

808 int z = i / (Y * X); 

809 int y = (i % (Y * X)) / X; 

810 int x = i % X; 

811 

812 // Calculate base index for this slice 

813 int slice_base = z * Y * X; 

814 

815 // Helper function to get pixel with configurable boundary handling 

816 auto get_pixel = [&](int py, int px) -> T { 

817 if (mode == 0) { // constant 

818 if (py < 0 || py >= Y || px < 0 || px >= X) return cval; 

819 } else if (mode == 1) { // reflect 

820 py = py < 0 ? -py : (py >= Y ? 2*Y - py - 2 : py); 

821 px = px < 0 ? -px : (px >= X ? 2*X - px - 2 : px); 

822 } else if (mode == 2) { // nearest 

823 py = py < 0 ? 0 : (py >= Y ? Y-1 : py); 

824 px = px < 0 ? 0 : (px >= X ? X-1 : px); 

825 } else if (mode == 3) { // wrap 

826 py = py < 0 ? py + Y : (py >= Y ? py - Y : py); 

827 px = px < 0 ? px + X : (px >= X ? px - X : px); 

828 } 

829 return input[slice_base + py * X + px]; 

830 }; 

831 

832 // Sobel X kernel: [[-1,0,1],[-2,0,2],[-1,0,1]] (within slice only) 

833 T gx = -get_pixel(y-1, x-1) + get_pixel(y-1, x+1) + 

834 -2*get_pixel(y, x-1) + 2*get_pixel(y, x+1) + 

835 -get_pixel(y+1, x-1) + get_pixel(y+1, x+1); 

836 

837 // Sobel Y kernel: [[-1,-2,-1],[0,0,0],[1,2,1]] (within slice only) 

838 T gy = -get_pixel(y-1, x-1) - 2*get_pixel(y-1, x) - get_pixel(y-1, x+1) + 

839 get_pixel(y+1, x-1) + 2*get_pixel(y+1, x) + get_pixel(y+1, x+1); 

840 

841 // Calculate magnitude 

842 output = sqrt(gx*gx + gy*gy); 

843 ''', 

844 'sobel_2d_parallel' 

845 ) 

846 return _sobel_2d_parallel 

847 

848@cupy_func 

849def sobel_2d_vectorized(image: "cp.ndarray", mode: str = "reflect", cval: float = 0.0) -> "cp.ndarray": 

850 """ 

851 Apply 2D Sobel edge detection to all slices simultaneously - TRUE GPU PARALLELIZED. 

852 

853 Each slice is treated as an independent 2D grayscale image. All pixels across 

854 all slices are processed in parallel on the GPU with slice independence guaranteed. 

855 Uses ElementwiseKernel for maximum performance and true parallelization. 

856 

857 Args: 

858 image: 3D CuPy array of shape (Z, Y, X) 

859 mode: Boundary handling mode ('constant', 'reflect', 'nearest', 'wrap') 

860 cval: Constant value for 'constant' mode 

861 

862 Returns: 

863 Edge magnitude as 3D CuPy array of shape (Z, Y, X) 

864 """ 

865 _validate_3d_array(image) 

866 

867 # Map mode string to integer 

868 mode_map = {'constant': 0, 'reflect': 1, 'nearest': 2, 'wrap': 3} 

869 if mode not in mode_map: 

870 raise ValueError(f"Unknown mode: {mode}. Valid modes: {list(mode_map.keys())}") 

871 

872 mode_int = mode_map[mode] 

873 Z, Y, X = image.shape 

874 input_float = image.astype(cp.float32) 

875 output = cp.zeros_like(input_float) 

876 

877 # Launch parallel kernel - each thread processes one pixel 

878 sobel_kernel = _get_sobel_2d_kernel() 

879 sobel_kernel(input_float, Y, X, mode_int, cp.float32(cval), output) 

880 

881 return output.astype(image.dtype) 

882 

883@cupy_func 

884def sobel_3d_voxel(image: "cp.ndarray", mode: str = "reflect", cval: float = 0.0) -> "cp.ndarray": 

885 """ 

886 Apply true 3D voxel Sobel edge detection including Z-axis gradients - GPU PARALLELIZED. 

887 

888 This computes gradients in all three spatial dimensions (X, Y, Z) for true 

889 volumetric edge detection. Useful for 3D structure analysis where Z-dimension 

890 has spatial meaning (not just independent slices). 

891 

892 Args: 

893 image: 3D CuPy array of shape (Z, Y, X) 

894 mode: Boundary handling mode ('constant', 'reflect', 'nearest', 'wrap', 'mirror') 

895 cval: Constant value for 'constant' mode 

896 

897 Returns: 

898 3D edge magnitude as 3D CuPy array of shape (Z, Y, X) 

899 """ 

900 _validate_3d_array(image) 

901 

902 # Convert to float32 for processing 

903 image_float = image.astype(cp.float32) 

904 

905 # Apply Sobel filters in all three directions - GPU PARALLELIZED 

906 sobel_x = ndimage.sobel(image_float, axis=2, mode=mode, cval=cval) # X-direction gradients 

907 sobel_y = ndimage.sobel(image_float, axis=1, mode=mode, cval=cval) # Y-direction gradients