|
| 1 | +"""Module to evaluate diagnostic metrics on tracks.""" |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +from scipy.spatial.distance import cdist |
| 5 | + |
| 6 | +from spine.ana.base import AnaBase |
| 7 | + |
| 8 | +from spine.utils.globals import TRACK_SHP |
| 9 | +from spine.utils.numba_local import principal_components |
| 10 | + |
| 11 | + |
| 12 | +__all__ = ['TrackCompletenessAna'] |
| 13 | + |
| 14 | + |
| 15 | +class TrackCompletenessAna(AnaBase): |
| 16 | + """This analysis script identifies gaps in tracks and measures the |
| 17 | + cumulative length of these gaps relative to the track length. |
| 18 | +
|
| 19 | + This is a useful diagnostic tool to evaluate the space-point efficiency |
| 20 | + on tracks (good standard candal as track should have exactly no gap in |
| 21 | + a perfectly efficient detector). |
| 22 | + """ |
| 23 | + |
| 24 | + # Name of the analysis script (as specified in the configuration) |
| 25 | + name = 'track_completeness' |
| 26 | + |
| 27 | + def __init__(self, time_window=None, run_mode='both', |
| 28 | + truth_point_mode='points', **kwargs): |
| 29 | + """Initialize the analysis script. |
| 30 | +
|
| 31 | + Parameters |
| 32 | + ---------- |
| 33 | + time_window : List[float] |
| 34 | + Time window within which to include particle (only works for `truth`) |
| 35 | + **kwargs : dict, optional |
| 36 | + Additional arguments to pass to :class:`AnaBase` |
| 37 | + """ |
| 38 | + # Initialize the parent class |
| 39 | + super().__init__('particle', run_mode, truth_point_mode, **kwargs) |
| 40 | + |
| 41 | + # Store the time window |
| 42 | + self.time_window = time_window |
| 43 | + assert time_window is None or len(time_window) == 2, ( |
| 44 | + "Time window must be specified as an array of two values.") |
| 45 | + assert time_window is None or run_mode == 'truth', ( |
| 46 | + "Time of reconstructed particle is unknown.") |
| 47 | + |
| 48 | + # Make sure the metadata is provided (rasterization needed) |
| 49 | + self.update_keys({'meta': True}) |
| 50 | + |
| 51 | + # Initialize the CSV writer(s) you want |
| 52 | + for prefix in self.prefixes: |
| 53 | + self.initialize_writer(prefix) |
| 54 | + |
| 55 | + def process(self, data): |
| 56 | + """Evaluate track completeness for tracks in one entry. |
| 57 | +
|
| 58 | + Parameters |
| 59 | + ---------- |
| 60 | + data : dict |
| 61 | + Dictionary of data products |
| 62 | + """ |
| 63 | + # Fetch the pixel size in this image (assume cubic cells) |
| 64 | + pixel_size = data['meta'].size[0] |
| 65 | + |
| 66 | + # Loop over the types of particle data products |
| 67 | + for key in self.obj_keys: |
| 68 | + # Fetch the prefix ('reco' or 'truth') |
| 69 | + prefix = key.split('_')[0] |
| 70 | + |
| 71 | + # Loop over particle objects |
| 72 | + for part in data[key]: |
| 73 | + # Check that the particle is a track |
| 74 | + if part.shape != TRACK_SHP: |
| 75 | + continue |
| 76 | + |
| 77 | + # If needed, check on the particle time |
| 78 | + if self.time_window is not None: |
| 79 | + if part.t < self.time_window[0] or part.t > self.time_window[1]: |
| 80 | + continue |
| 81 | + |
| 82 | + # Initialize the particle dictionary |
| 83 | + comp_dict = {'particle_id': part.id} |
| 84 | + |
| 85 | + # Fetch the particle point coordinates |
| 86 | + points = self.get_points(part) |
| 87 | + |
| 88 | + # Find start/end points, collapse onto track cluster |
| 89 | + start = points[np.argmin(cdist([part.start_point], points))] |
| 90 | + end = points[np.argmin(cdist([part.end_point], points))] |
| 91 | + |
| 92 | + # Add the direction of the track |
| 93 | + vec = end - start |
| 94 | + length = np.linalg.norm(vec) |
| 95 | + if length: |
| 96 | + vec /= length |
| 97 | + |
| 98 | + comp_dict['size'] = len(points) |
| 99 | + comp_dict['length'] = length |
| 100 | + comp_dict.update( |
| 101 | + {'dir_x': vec[0], 'dir_y': vec[1], 'dir_z': vec[2]}) |
| 102 | + |
| 103 | + # Chunk out the track along gaps, estimate gap length |
| 104 | + chunk_labels = self.cluster_track_chunks( |
| 105 | + points, start, end, pixel_size) |
| 106 | + gaps = self.sequential_cluster_distances( |
| 107 | + points, chunk_labels, start) |
| 108 | + |
| 109 | + # Substract minimum gap distance due to rasterization |
| 110 | + min_gap = pixel_size/np.max(np.abs(vec)) |
| 111 | + gaps -= min_gap |
| 112 | + |
| 113 | + # Store gap information |
| 114 | + comp_dict['num_gaps'] = len(gaps) |
| 115 | + comp_dict['gap_length'] = np.sum(gaps) |
| 116 | + comp_dict['gap_frac'] = np.sum(gaps)/length |
| 117 | + |
| 118 | + # Append the dictionary to the CSV log |
| 119 | + self.append(prefix, **comp_dict) |
| 120 | + |
| 121 | + @staticmethod |
| 122 | + def cluster_track_chunks(points, start_point, end_point, pixel_size): |
| 123 | + """Find point where the track is broken, divide out the track |
| 124 | + into self-contained chunks which are Linf connect (Moore neighbors). |
| 125 | +
|
| 126 | + Parameters |
| 127 | + ---------- |
| 128 | + points : np.ndarray |
| 129 | + (N, 3) List of track cluster point coordinates |
| 130 | + start_point : np.ndarray |
| 131 | + (3) Start point of the track cluster |
| 132 | + end_point : np.ndarray |
| 133 | + (3) End point of the track cluster |
| 134 | + pixel_size : float |
| 135 | + Dimension of one pixel, used to identify what is big enough to |
| 136 | + constitute a break |
| 137 | +
|
| 138 | + Returns |
| 139 | + ------- |
| 140 | + np.ndarray |
| 141 | + (N) Track chunk labels |
| 142 | + """ |
| 143 | + # Project and cluster on the projected axis |
| 144 | + direction = (end_point-start_point)/np.linalg.norm(end_point-start_point) |
| 145 | + projs = np.dot(points - start_point, direction) |
| 146 | + perm = np.argsort(projs) |
| 147 | + seps = projs[perm][1:] - projs[perm][:-1] |
| 148 | + breaks = np.where(seps > pixel_size*1.1)[0] + 1 |
| 149 | + cluster_labels = np.empty(len(projs), dtype=int) |
| 150 | + for i, index in enumerate(np.split(np.arange(len(projs)), breaks)): |
| 151 | + cluster_labels[perm[index]] = i |
| 152 | + |
| 153 | + return cluster_labels |
| 154 | + |
| 155 | + @staticmethod |
| 156 | + def sequential_cluster_distances(points, labels, start_point): |
| 157 | + """Order clusters in order of distance from a starting point, compute |
| 158 | + the distances between successive clusters. |
| 159 | +
|
| 160 | + Parameters |
| 161 | + ---------- |
| 162 | + points : np.ndarray |
| 163 | + (N, 3) List of track cluster point coordinates |
| 164 | + labels : np.ndarray |
| 165 | + (N) Track chunk labels |
| 166 | + start_point : np.ndarray |
| 167 | + (3) Start point of the track cluster |
| 168 | + """ |
| 169 | + # If there's only one cluster, nothing to do here |
| 170 | + unique_labels = np.unique(labels) |
| 171 | + if len(unique_labels) < 2: |
| 172 | + return np.empty(0, dtype=float), np.empty(0, dtype=float) |
| 173 | + |
| 174 | + # Order clusters |
| 175 | + start_dist = cdist([start_point], points).flatten() |
| 176 | + start_clust_dist = np.empty(len(unique_labels)) |
| 177 | + for i, c in enumerate(unique_labels): |
| 178 | + start_clust_dist[i] = np.min(start_dist[labels == c]) |
| 179 | + ordered_labels = unique_labels[np.argsort(start_clust_dist)] |
| 180 | + |
| 181 | + # Compute the intercluster distance and relative angle |
| 182 | + n_gaps = len(ordered_labels) - 1 |
| 183 | + dists = np.empty(n_gaps, dtype=float) |
| 184 | + for i in range(n_gaps): |
| 185 | + points_i = points[labels == ordered_labels[i]] |
| 186 | + points_j = points[labels == ordered_labels[i + 1]] |
| 187 | + dists[i] = np.min(cdist(points_i, points_j)) |
| 188 | + |
| 189 | + return dists |
0 commit comments