6279221893

Committed 22 Sep 2023 09:19PM UTC coverage: 74.349% (+0.2%) from 74.109%

Build # 6279221893

Build Type

Pull #19

github

Committed by

WardLT

Commit Message

Keep track whether voids touch sides in tracks

Pull Request Pull Request #19: Keep track of whether voids touch sides in tracks

Run Details

7 of 7 new or added lines in 1 file covered. (100.0%)

400 of 538 relevant lines covered (74.35%)

0.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.56

/rtdefects/analysis.py

"""Functions to analyze segmented images"""
import logging
from typing import Iterator, Tuple

from skimage import measure, morphology
from scipy.stats import siegelslopes
from scipy.interpolate import interp1d
import pandas as pd
import numpy as np

logger = logging.getLogger(__name__)


def analyze_defects(mask: np.ndarray, min_size: int = 50, edge_buffer: int = 8) -> Tuple[dict, np.ndarray]:
    """Analyze the voids in a masked image

    Args:
        mask: Masks for a defect image
        min_size: Minimum size of defects
        edge_buffer: Label voids as touching edge if they are within this many pixels of the side
    Returns:
        - Dictionary of the computed properties
        - Labeled images
    """

    # Clean up the mask
    mask = morphology.remove_small_objects(mask, min_size=min_size)
    mask = morphology.remove_small_holes(mask, min_size)
    mask = morphology.binary_erosion(mask, morphology.square(1))
    output = {'void_frac': mask.sum() / (mask.shape[0] * mask.shape[1])}

    # Assign labels to the labeled regions
    labels = measure.label(mask)
    output['void_count'] = int(labels.max())

    # Compute region properties
    props = measure.regionprops(labels, mask)
    radii = [p['equivalent_diameter'] / 2 for p in props]
    output['radii'] = radii
    output['radii_average'] = np.average(radii)
    output['positions'] = [p['centroid'][::-1] for p in props]  # From (y, x) to (x, y)

    # Determine if it touches the side
    output['touches_side'] = [
        min(p['bbox']) <= edge_buffer
        or p['bbox'][2] >= mask.shape[0] - edge_buffer
        or p['bbox'][3] >= mask.shape[1] - edge_buffer
        for p in props
    ]

    return output, labels


def convert_to_per_particle(per_frame: pd.DataFrame, position_col: str = 'positions') -> Iterator[pd.DataFrame]:
    """Convert the per-frame void information to the per-particle format expected by trackpy

    Args:
        per_frame: A DataFrame where each row is a different image and
            contains the defect locations in `positions` and sizes in `radii` columns.
        position_col: Name of the column holding positions of the particles
    Yields:
        A dataframe where each row is a different defect
    """

    for rid, row in per_frame.iterrows():
        particles = pd.DataFrame(row[position_col], columns=['x', 'y'])
        particles['local_id'] = np.arange(len(row['positions']))
        particles['frame'] = rid
        particles['radius'] = row['radii']
        particles['touches_side'] = row['touches_side']
        yield particles


def compute_drift(tracks: pd.DataFrame, minimum_tracks: int = 1) -> np.ndarray:
    """Estimate the drift for each frame from the positions of voids that were mapped between multiple frames

    We determine the "drift" based on the median displacement of all voids, which is based
    on the assumption that there is no net motion of all the voids.

    We compute the drift in each frame and assume the drift remains unchanged if there are no voids matched
    between a frame and the previous.

    In contrast, trackpy uses the mean and only computes drift when there are matches between frames.

    Args:
        tracks: Track information generated by trackpy.
        minimum_tracks: The minimum number of tracks a void must appear to be used in drift correction
    Returns:
        Drift correction for each frame
    """

    # We'll assume that the first frame has a void
    drifts = [(0, 0)]

    # We're going to go frame-by-frame and guess the drift from the previous frame
    last_frame = tracks.query('frame==0')
    for fid in range(1, tracks['frame'].max() + 1):
        # Join the two frames
        my_frame = tracks.query(f'frame=={fid}')
        aligned = last_frame.merge(my_frame, on='particle')

        # The current frame will be the previous for the next iteration
        last_frame = my_frame

        # If there are no voids in both frames, assign a drift change of 0
        if len(aligned) < minimum_tracks:
            drifts.append(drifts[-1])
            continue

        # Get the median displacements displacements
        last_pos = aligned[['x_x', 'y_x']].values
        cur_pos = aligned[['x_y', 'y_y']].values
        median_disp = np.mean(cur_pos - last_pos, axis=0)

        # Add the drift to that of the previous image
        drift = np.add(drifts[-1], median_disp)
        drifts.append(drift)

    return np.array(drifts)


def compile_void_tracks(tracks: pd.DataFrame) -> pd.DataFrame:
    """Compile summary statistics about each void

    Args:
        tracks: Track information for each void over time

    Returns:
        Dataframe of the summary of voids
        - "start_frame": First frame in which the void appears
        - "end_frame": Last frame in which the void appears
        - "total_frames": Total number of frames in which the void appears
        - "positions": Positions of the void in each frame
        - "touches_side": Whether the void touches the side at this frame
        - "local_id": ID of the void in each frame (if available)
        - "disp_from_start": How far the void has moved from the first frame
        - "max_disp": Maximum distance the void moved
        - "drift_rate": Average displacement from center over time
        - "dist_traveled": Total path distance the void has traveled
        - "total_travel": How far the void traveled over its whole life
        - "movement_rate": How far the void moves per frame
        - "radii": Radius of the void in each frame
        - "max_radius": Maximum radius of the void
        - "min_radius": Minimum radius of the void
        - "growth_rate": Median rate of change of the radius
    """

    # Loop over all unique voids
    voids = []
    for t, track in tracks.groupby('particle'):
        # Get the frames where this void is visible
        visible_frames = track['frame']

        # Get all frames between start and stop
        frames_id = np.arange(track['frame'].min(), track['frame'].max() + 1)

        # Build an interpolator for position as a function of frame
        if len(track) == 1:
            positions = track[['x', 'y']].values
        else:
            x_inter = interp1d(track['frame'], track['x'])
            y_inter = interp1d(track['frame'], track['y'])

            # Compute the displacement over each step
            positions = [(x_inter(f), y_inter(f)) for f in frames_id]
            positions = np.array(positions)

        # Get the ID from each frame and whether it touches the side
        id_lookup = dict(zip(track['frame'], track['local_id']))
        local_id = [id_lookup.get(i, None) for i in frames_id]

        # Use interpolation to detect if any point used to interpolate
        #  the void position was on the side
        if len(track) == 1:
            touches_side = track['touches_side'].values
        else:
            ts_inter = interp1d(track['frame'], np.array(track['touches_side'], dtype=float), kind='linear')
            touches_side = ts_inter(frames_id) > 0  # It is only zero if neither point used in the left or right touches side (and equals 1)

        # Gather some basic information about the void
        void_info = {
            'start_frame': np.min(visible_frames),
            'end_frame': np.max(visible_frames),
            'total_frames': len(frames_id),
            'inferred_frames': len(frames_id) - len(track),
            'positions': positions,
            'touches_side': touches_side,
            'local_id': local_id
        }

        # If there is only one frame, we cannot do the following steps
        if positions.shape[0] > 1:
            # Compute the displacement from the start
            void_info['disp_from_start'] = np.linalg.norm(positions - positions[0, :], axis=1)
            void_info['max_disp'] = np.max(void_info['disp_from_start'])
            void_info['drift_rate'] = void_info['max_disp'] / void_info['total_frames']

            # Get the displacement for each step
            void_info['dist_traveled'] = np.concatenate(([0], np.cumsum(np.linalg.norm(np.diff(positions, axis=0), axis=1))))
            void_info['total_traveled'] = void_info['dist_traveled'][-1]
            void_info['movement_rate'] = void_info['total_traveled'] / void_info['total_frames']

        # More stats if we have radii
        if len(track) == 1:
            radii = track['radius'].values
        else:
            r_inter = interp1d(track['frame'], track['radius'])
            radii = r_inter(frames_id)

        # Store some summary information
        void_info['radii'] = radii
        void_info['max_radius'] = max(radii)
        void_info['min_radius'] = min(radii)
        if len(radii) > 3:
            void_info['growth_rate'] = siegelslopes(radii)[0]

        # Add it to list
        voids.append(void_info)
    return pd.DataFrame(voids)

1	"""Functions to analyze segmented images"""
2	import logging	1✔
3	from typing import Iterator, Tuple	1✔
4
5	from skimage import measure, morphology	1✔
6	from scipy.stats import siegelslopes	1✔
7	from scipy.interpolate import interp1d	1✔
8	import pandas as pd	1✔
9	import numpy as np	1✔
10
11	logger = logging.getLogger(__name__)	1✔
12
13
14	def analyze_defects(mask: np.ndarray, min_size: int = 50, edge_buffer: int = 8) -> Tuple[dict, np.ndarray]:	1✔
15	"""Analyze the voids in a masked image
16
17	Args:
18	mask: Masks for a defect image
19	min_size: Minimum size of defects
20	edge_buffer: Label voids as touching edge if they are within this many pixels of the side
21	Returns:
22	- Dictionary of the computed properties
23	- Labeled images
24	"""
25
26	# Clean up the mask
27	mask = morphology.remove_small_objects(mask, min_size=min_size)	1✔
28	mask = morphology.remove_small_holes(mask, min_size)	1✔
29	mask = morphology.binary_erosion(mask, morphology.square(1))	1✔
30	output = {'void_frac': mask.sum() / (mask.shape[0] * mask.shape[1])}	1✔
31
32	# Assign labels to the labeled regions
33	labels = measure.label(mask)	1✔
34	output['void_count'] = int(labels.max())	1✔
35
36	# Compute region properties
37	props = measure.regionprops(labels, mask)	1✔
38	radii = [p['equivalent_diameter'] / 2 for p in props]	1✔
39	output['radii'] = radii	1✔
40	output['radii_average'] = np.average(radii)	1✔
41	output['positions'] = [p['centroid'][::-1] for p in props] # From (y, x) to (x, y)	1✔
42
43	# Determine if it touches the side
44	output['touches_side'] = [	1✔
45	min(p['bbox']) <= edge_buffer
46	or p['bbox'][2] >= mask.shape[0] - edge_buffer
47	or p['bbox'][3] >= mask.shape[1] - edge_buffer
48	for p in props
49	]
50
51	return output, labels	1✔
52
53
54	def convert_to_per_particle(per_frame: pd.DataFrame, position_col: str = 'positions') -> Iterator[pd.DataFrame]:	1✔
55	"""Convert the per-frame void information to the per-particle format expected by trackpy
56
57	Args:
58	per_frame: A DataFrame where each row is a different image and
59	contains the defect locations in `positions` and sizes in `radii` columns.
60	position_col: Name of the column holding positions of the particles
61	Yields:
62	A dataframe where each row is a different defect
63	"""
64
65	for rid, row in per_frame.iterrows():	1✔
66	particles = pd.DataFrame(row[position_col], columns=['x', 'y'])	1✔
67	particles['local_id'] = np.arange(len(row['positions']))	1✔
68	particles['frame'] = rid	1✔
69	particles['radius'] = row['radii']	1✔
70	particles['touches_side'] = row['touches_side']	1✔
71	yield particles	1✔
72
73
74	def compute_drift(tracks: pd.DataFrame, minimum_tracks: int = 1) -> np.ndarray:	1✔
75	"""Estimate the drift for each frame from the positions of voids that were mapped between multiple frames
76
77	We determine the "drift" based on the median displacement of all voids, which is based
78	on the assumption that there is no net motion of all the voids.
79
80	We compute the drift in each frame and assume the drift remains unchanged if there are no voids matched
81	between a frame and the previous.
82
83	In contrast, trackpy uses the mean and only computes drift when there are matches between frames.
84
85	Args:
86	tracks: Track information generated by trackpy.
87	minimum_tracks: The minimum number of tracks a void must appear to be used in drift correction
88	Returns:
89	Drift correction for each frame
90	"""
91
92	# We'll assume that the first frame has a void
93	drifts = [(0, 0)]	1✔
94
95	# We're going to go frame-by-frame and guess the drift from the previous frame
96	last_frame = tracks.query('frame==0')	1✔
97	for fid in range(1, tracks['frame'].max() + 1):	1✔
98	# Join the two frames
99	my_frame = tracks.query(f'frame=={fid}')	1✔
100	aligned = last_frame.merge(my_frame, on='particle')	1✔
101
102	# The current frame will be the previous for the next iteration
103	last_frame = my_frame	1✔
104
105	# If there are no voids in both frames, assign a drift change of 0
106	if len(aligned) < minimum_tracks:	1✔
107	drifts.append(drifts[-1])	×
108	continue	×
109
110	# Get the median displacements displacements
111	last_pos = aligned[['x_x', 'y_x']].values	1✔
112	cur_pos = aligned[['x_y', 'y_y']].values	1✔
113	median_disp = np.mean(cur_pos - last_pos, axis=0)	1✔
114
115	# Add the drift to that of the previous image
116	drift = np.add(drifts[-1], median_disp)	1✔
117	drifts.append(drift)	1✔
118
119	return np.array(drifts)	1✔
120
121
122	def compile_void_tracks(tracks: pd.DataFrame) -> pd.DataFrame:	1✔
123	"""Compile summary statistics about each void
124
125	Args:
126	tracks: Track information for each void over time
127
128	Returns:
129	Dataframe of the summary of voids
130	- "start_frame": First frame in which the void appears
131	- "end_frame": Last frame in which the void appears
132	- "total_frames": Total number of frames in which the void appears
133	- "positions": Positions of the void in each frame
134	- "touches_side": Whether the void touches the side at this frame
135	- "local_id": ID of the void in each frame (if available)
136	- "disp_from_start": How far the void has moved from the first frame
137	- "max_disp": Maximum distance the void moved
138	- "drift_rate": Average displacement from center over time
139	- "dist_traveled": Total path distance the void has traveled
140	- "total_travel": How far the void traveled over its whole life
141	- "movement_rate": How far the void moves per frame
142	- "radii": Radius of the void in each frame
143	- "max_radius": Maximum radius of the void
144	- "min_radius": Minimum radius of the void
145	- "growth_rate": Median rate of change of the radius
146	"""
147
148	# Loop over all unique voids
149	voids = []	1✔
150	for t, track in tracks.groupby('particle'):	1✔
151	# Get the frames where this void is visible
152	visible_frames = track['frame']	1✔
153
154	# Get all frames between start and stop
155	frames_id = np.arange(track['frame'].min(), track['frame'].max() + 1)	1✔
156
157	# Build an interpolator for position as a function of frame
158	if len(track) == 1:	1✔
159	positions = track[['x', 'y']].values	1✔
160	else:
161	x_inter = interp1d(track['frame'], track['x'])	1✔
162	y_inter = interp1d(track['frame'], track['y'])	1✔
163
164	# Compute the displacement over each step
165	positions = [(x_inter(f), y_inter(f)) for f in frames_id]	1✔
166	positions = np.array(positions)	1✔
167
168	# Get the ID from each frame and whether it touches the side
169	id_lookup = dict(zip(track['frame'], track['local_id']))	1✔
170	local_id = [id_lookup.get(i, None) for i in frames_id]	1✔
171
172	# Use interpolation to detect if any point used to interpolate
173	# the void position was on the side
174	if len(track) == 1:	1✔
175	touches_side = track['touches_side'].values	1✔
176	else:
177	ts_inter = interp1d(track['frame'], np.array(track['touches_side'], dtype=float), kind='linear')	1✔
178	touches_side = ts_inter(frames_id) > 0 # It is only zero if neither point used in the left or right touches side (and equals 1)	1✔
179
180	# Gather some basic information about the void
181	void_info = {	1✔
182	'start_frame': np.min(visible_frames),
183	'end_frame': np.max(visible_frames),
184	'total_frames': len(frames_id),
185	'inferred_frames': len(frames_id) - len(track),
186	'positions': positions,
187	'touches_side': touches_side,
188	'local_id': local_id
189	}
190
191	# If there is only one frame, we cannot do the following steps
192	if positions.shape[0] > 1:	1✔
193	# Compute the displacement from the start
194	void_info['disp_from_start'] = np.linalg.norm(positions - positions[0, :], axis=1)	1✔
195	void_info['max_disp'] = np.max(void_info['disp_from_start'])	1✔
196	void_info['drift_rate'] = void_info['max_disp'] / void_info['total_frames']	1✔
197
198	# Get the displacement for each step
199	void_info['dist_traveled'] = np.concatenate(([0], np.cumsum(np.linalg.norm(np.diff(positions, axis=0), axis=1))))	1✔
200	void_info['total_traveled'] = void_info['dist_traveled'][-1]	1✔
201	void_info['movement_rate'] = void_info['total_traveled'] / void_info['total_frames']	1✔
202
203	# More stats if we have radii
204	if len(track) == 1:	1✔
205	radii = track['radius'].values	1✔
206	else:
207	r_inter = interp1d(track['frame'], track['radius'])	1✔
208	radii = r_inter(frames_id)	1✔
209
210	# Store some summary information
211	void_info['radii'] = radii	1✔
212	void_info['max_radius'] = max(radii)	1✔
213	void_info['min_radius'] = min(radii)	1✔
214	if len(radii) > 3:	1✔
215	void_info['growth_rate'] = siegelslopes(radii)[0]	1✔
216
217	# Add it to list
218	voids.append(void_info)	1✔
219	return pd.DataFrame(voids)	1✔

ivem-argonne / real-time-defect-analysis / 6279221893

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous