Coverage for src/pyhiperta/cleaning.py: 100%

1# Copyright 2026 CNRS 

2# This software is distributed under the terms of the CeCILL-C free software license. 

3 

4"""Charge images cleaning algorithms: tail-cut cleaning.""" 

5 

6import numpy as np 

7 

8from pyhiperta.utils.convolve import convolve_view 

9from pyhiperta.waveform_indexing import neighbors_only_stencil 

10 

11 

12def tail_cuts_cleaning( 

13 waveforms_2D: np.ndarray, 

14 pixel_threshold: float, 

15 neighbors_threshold: float, 

16 min_number_neighbors: int, 

17) -> np.ndarray: 

18 """Compute the mask of "signal" pixel that pass the tail_cuts thresholds. 

19 

20 The implementation is in 2 steps: 

21 - find the group of pixels that pass the `pixel_threshold` 

22 - find the pixels that pass `neighbors_threshold` and have at least 1 neighbor passing the 1st step. 

23 

24 Parameters 

25 ---------- 

26 waveforms_2D : np.ndarray 

27 Batch or integrated waveform in 2D format. Shape: ([N_batch,] N_pixels_x, N_pixels_y) 

28 pixel_threshold : float 

29 A pixel with a value greater or equal than `pixel_threshold` and at least `min_number_neighbors` neighbors 

30 that have a value greater or equal than `pixel_threshold` are considered "signal". 

31 neighbors_threshold : float 

32 A pixel with a value greater or equal than `neighbors_threshold` and at least 1 neighbor that is considered 

33 signal according to `pixel_threshold` will be considered "signal" as well. 

34 min_number_neighbors : int 

35 Minimum number of neighboring pixels that must have a value above `pixel_threshold` to be considered "signal". 

36 

37 Returns 

38 ------- 

39 np.ndarray 

40 A boolean mask with value True for "signal" pixels and value False otherwise. 

41 

42 Raises 

43 ------ 

44 ValueError 

45 If the shape of waveforms_2D can not be interpreted as a (batch of) 2D waveforms 

46 """ 

47 if len(waveforms_2D.shape) < 2: 

48 raise ValueError( 

49 "waveforms must be an array with at least 2 dimensions " 

50 f"(waveform 2D or batch of waveform 2D), but got {waveforms_2D.shape}" 

51 ) 

52 

53 # kepp pixels that are 

54 # 1: above pixel threshold and have at least min_number_neighbors above pixel threshold as well 

55 # 2: pixels that are above neighbor's threshold and have at least 1 neighbor that checks condition 1 

56 

57 nb_batch_dimension = len(waveforms_2D.shape) - 2 

58 

59 neighbors_stencil = neighbors_only_stencil() 

60 # add as many dimension to the 2D neighbor stencil as required (to allow for batch dimension) 

61 neighbors_stencil = neighbors_stencil[*([np.newaxis] * nb_batch_dimension), ...] 

62 # get the axis dimension to reduce when reducing the view: 

63 # If waveform2D.shape = (3, 55, 55) then stencil will have shape (1, 3, 3) and the 

64 # view will have shape (3, 55, 55, 1, 3, 3) 

65 # To reduce the view (compute the convolved operation) we will reduce on axis -3, -2, -1 

66 convolution_reduction_axis = tuple([-i - 1 for i in range(len(neighbors_stencil.shape))]) 

67 

68 # we will pad with one 0 on both ends of waveform 2D, and not pad the remaining (batch) axis 

69 pad_values = [(0, 0)] * nb_batch_dimension + [(1, 1), (1, 1)] 

70 

71 # Pad the waveform with 0 on all edges to be able to convolve without reducing the shape 

72 waveforms_2D_padded = np.pad(waveforms_2D, pad_values, mode="constant", constant_values=0) 

73 neighbors_only_view = convolve_view(waveforms_2D_padded, neighbors_stencil.shape) * neighbors_stencil 

74 # condition 1: 

75 mask = (waveforms_2D >= pixel_threshold) & ( 

76 (neighbors_only_view >= pixel_threshold).sum(axis=convolution_reduction_axis) >= min_number_neighbors 

77 ) 

78 # pad the mask to compute condition 2: condition on neighbor's number of neighbors 

79 padded_mask = np.pad(mask, pad_values, mode="constant", constant_values=0) 

80 # get the convolution view for the neighbors passing condition 1 

81 neighbors_passing_condition_1 = convolve_view(padded_mask, neighbors_stencil.shape) * neighbors_stencil 

82 # condition 2: 

83 mask |= (waveforms_2D >= neighbors_threshold) & ( 

84 neighbors_passing_condition_1.any(axis=convolution_reduction_axis) 

85 ) 

86 return mask