Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Directly return peak position from neighbor average waveform #1914

Merged
merged 1 commit into from
May 24, 2022
Merged

Directly return peak position from neighbor average waveform #1914

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions ctapipe/image/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@
TwoPassWindowSum,
extract_around_peak,
extract_sliding_window,
neighbor_average_waveform,
neighbor_average_maximum,
subtract_baseline,
integration_correction,
)
Expand Down Expand Up @@ -110,7 +110,7 @@
"TwoPassWindowSum",
"extract_around_peak",
"extract_sliding_window",
"neighbor_average_waveform",
"neighbor_average_maximum",
"subtract_baseline",
"integration_correction",
"DataVolumeReducer",
Expand Down
18 changes: 8 additions & 10 deletions ctapipe/image/extractor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
"TwoPassWindowSum",
"extract_around_peak",
"extract_sliding_window",
"neighbor_average_waveform",
"neighbor_average_maximum",
"subtract_baseline",
"integration_correction",
]
Expand Down Expand Up @@ -195,7 +195,7 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time)


@njit(cache=True)
def neighbor_average_waveform(waveforms, neighbors_indices, neighbors_indptr, lwt):
def neighbor_average_maximum(waveforms, neighbors_indices, neighbors_indptr, lwt):
"""
Obtain the average waveform built from the neighbors of each pixel

Expand Down Expand Up @@ -228,19 +228,18 @@ def neighbor_average_waveform(waveforms, neighbors_indices, neighbors_indptr, lw

# initialize to waveforms weighted with lwt
# so the value of the pixel itself is already taken into account
average = waveforms * lwt
peak_pos = np.empty(n_pixels, dtype=np.int64)

for pixel in prange(n_pixels):
average = waveforms[pixel] * lwt
neighbors = indices[indptr[pixel] : indptr[pixel + 1]]

n = lwt
for neighbor in neighbors:
average[pixel] += waveforms[neighbor]
n += 1
average += waveforms[neighbor]

average[pixel] /= n
peak_pos[pixel] = np.argmax(average)

return average
return peak_pos


def subtract_baseline(waveforms, baseline_start, baseline_end):
Expand Down Expand Up @@ -751,13 +750,12 @@ def _calculate_correction(self, telid):

def __call__(self, waveforms, telid, selected_gain_channel):
neighbors = self.subarray.tel[telid].camera.geometry.neighbor_matrix_sparse
average_wfs = neighbor_average_waveform(
peak_index = neighbor_average_maximum(
waveforms,
neighbors_indices=neighbors.indices,
neighbors_indptr=neighbors.indptr,
lwt=self.lwt.tel[telid],
)
peak_index = average_wfs.argmax(axis=-1)
charge, peak_time = extract_around_peak(
waveforms,
peak_index,
Expand Down
16 changes: 9 additions & 7 deletions ctapipe/image/tests/test_extractor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
extract_around_peak,
extract_sliding_window,
integration_correction,
neighbor_average_waveform,
neighbor_average_maximum,
subtract_baseline,
)
from ctapipe.image.toymodel import SkewedGaussian, WaveformModel, obtain_time_image
Expand Down Expand Up @@ -191,10 +191,10 @@ def test_extract_around_peak_charge_expected(toymodel):
assert_equal(charge, n_samples)


def test_neighbor_average_waveform(toymodel):
def test_neighbor_average_peakpos(toymodel):
waveforms, subarray, telid, _, _, _ = toymodel
neighbors = subarray.tel[telid].camera.geometry.neighbor_matrix_sparse
average_wf = neighbor_average_waveform(
peak_pos = neighbor_average_maximum(
waveforms,
neighbors_indices=neighbors.indices,
neighbors_indptr=neighbors.indptr,
Expand All @@ -204,10 +204,11 @@ def test_neighbor_average_waveform(toymodel):
pixel = 0
_, nei_pixel = np.where(neighbors[pixel].A)
expected_average = waveforms[nei_pixel].sum(0) / len(nei_pixel)
assert_allclose(average_wf[pixel], expected_average, rtol=1e-3)
expected_peak_pos = np.argmax(expected_average, axis=-1)
assert (peak_pos[pixel] == expected_peak_pos).all()

lwt = 4
average_wf = neighbor_average_waveform(
peak_pos = neighbor_average_maximum(
waveforms,
neighbors_indices=neighbors.indices,
neighbors_indptr=neighbors.indptr,
Expand All @@ -218,7 +219,8 @@ def test_neighbor_average_waveform(toymodel):
_, nei_pixel = np.where(neighbors[pixel].A)
nei_pixel = np.concatenate([nei_pixel, [pixel] * lwt])
expected_average = waveforms[nei_pixel].sum(0) / len(nei_pixel)
assert_allclose(average_wf[pixel], expected_average, rtol=1e-3)
expected_peak_pos = np.argmax(expected_average, axis=-1)
assert (peak_pos[pixel] == expected_peak_pos).all()


def test_extract_peak_time_within_range():
Expand Down Expand Up @@ -532,7 +534,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray):

tel_id = 1
camera = subarray.tel[tel_id].camera
sample_rate = camera.readout.sampling_rate.to_value(u.ns ** -1)
sample_rate = camera.readout.sampling_rate.to_value(u.ns**-1)
n_pixels = camera.geometry.n_pixels
selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8)

Expand Down