Skip to content

Commit

Permalink
Merge pull request #1698 from HealthyPear/tests
Browse files Browse the repository at this point in the history 
Fix bugs in TwoPassWindowSum and define better unit-test
  • Loading branch information
@maxnoe
maxnoe authored May 10, 2021
2 parents 3af1488 + f8f4fdb commit b4d1e29
Show file tree
Hide file tree
Showing 2 changed files with 195 additions and 114 deletions.
163 changes: 82 additions & 81 deletions ctapipe/image/extractor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
from scipy.ndimage.filters import convolve1d
from typing import Tuple

from . import number_of_islands, largest_island, tailcuts_clean
from . import number_of_islands, tailcuts_clean, brightest_island
from .timing import timing_parameters
from .hillas import hillas_parameters, camera_to_shower_coordinates

Expand Down Expand Up @@ -781,7 +781,7 @@ class TwoPassWindowSum(ImageExtractor):
No information from neighboouring pixels is used.
#. Preliminary image cleaning via simple tailcut with minimum number
of core neighbours set at 1,
#. Only the biggest cluster of pixels is kept.
#. Only the brightest cluster of pixels is kept.
#. Parametrize following Hillas approach only if the resulting image has 3
or more pixels.
#. Do a linear fit of pulse time vs. distance along major image axis
Expand Down Expand Up @@ -994,50 +994,41 @@ def _apply_second_pass(

# STEP 3

# Indexes of pixels that will need the 2nd pass
non_core_pixels_mask = charge_1stpass < core_th

# find all islands using this cleaning
num_islands, labels = number_of_islands(camera_geometry, mask_clean)

if num_islands > 0:
# ...find the biggest one
mask_biggest = largest_island(labels)
# ...find the brightest one
mask_brightest_island = brightest_island(num_islands, labels, charge_1stpass)
else:
mask_biggest = mask_clean
mask_brightest_island = mask_clean

# for all pixels except the core ones in the main island of the
# preliminary image, the waveform will be integrated once more (2nd pass)

mask_2nd_pass = ~mask_brightest_island | (mask_brightest_island & (charge_1stpass < core_th))

# STEP 4

# if the resulting image has less then 3 pixels
# or there are more than 3 pixels but all contain a number of
# photoelectrons above the core threshold
if np.count_nonzero(mask_biggest) < 3:
if np.count_nonzero(mask_brightest_island) < 3:
# we return the 1st pass information
# NOTE: In this case, the image was not bright enough!
# We should label it as "bad and NOT use it"
return charge_1stpass, pulse_time_1stpass
elif np.count_nonzero(non_core_pixels_mask) == 0:
# Since all reconstructed charges are above the core threshold,
# there is no need to perform the 2nd pass.
# We return the 1st pass information.
# NOTE: In this case, even if this is 1st pass information,
# the image is actually very bright! We should label it as "good"!
return charge_1stpass, pulse_time_1stpass

# otherwise we proceed by parametrizing the image
camera_geometry_biggest = camera_geometry[mask_biggest]
charge_biggest = charge_1stpass[mask_biggest]
hillas = hillas_parameters(camera_geometry_biggest, charge_biggest)
camera_geometry_brightest = camera_geometry[mask_brightest_island]
charge_brightest = charge_1stpass[mask_brightest_island]
hillas = hillas_parameters(camera_geometry_brightest, charge_brightest)

# STEP 5

# linear fit of pulse time vs. distance along major image axis
# using only the main island surviving the preliminary
# image cleaning
timing = timing_parameters(
geom=camera_geometry_biggest,
image=charge_biggest,
peak_time=pulse_time_1stpass[mask_biggest],
geom=camera_geometry_brightest,
image=charge_brightest,
peak_time=pulse_time_1stpass[mask_brightest_island],
hillas_parameters=hillas,
)

Expand All @@ -1052,7 +1043,7 @@ def _apply_second_pass(

# get the predicted times as a linear relation
predicted_pulse_times = (
timing.slope * longitude[non_core_pixels_mask] + timing.intercept
timing.slope * longitude[mask_2nd_pass] + timing.intercept
)

# Convert time in ns to sample index using the sampling rate from
Expand All @@ -1064,14 +1055,13 @@ def _apply_second_pass(
predicted_pulse_times.value * self.sampling_rate_ghz[telid]
).astype(np.int64)

# Due to the fit these peak indexes can now be also outside of the
# readout window, so later we check for this.
# Due to the fit these peak indexes can lead to an integration window
# outside the readout window, so in the next step we check for this.

# STEP 6

# select only the waveforms correspondent to the non-core pixels
# of the main island survived from the 1st pass image cleaning
non_core_waveforms = waveforms[non_core_pixels_mask]
# select all pixels except the core ones in the main island
waveforms_to_repass = waveforms[mask_2nd_pass]

# Build 'width' and 'shift' arrays that adapt on the position of the
# window along each waveform
Expand All @@ -1080,44 +1070,58 @@ def _apply_second_pass(
# on the peak
window_width_default = 5
window_shift_default = 2

# now let's deal with some edge cases: the predicted peak falls before
# or after the readout window:
peak_before_readout = predicted_peaks < 0
peak_after_readout = predicted_peaks > (non_core_waveforms.shape[1] - 1)

# BUT, if the resulting 5-samples window falls outside of the readout
# window then we take the first (or last) 5 samples

# first we find where the integration window edges WOULD BE
integration_windows_start = predicted_peaks - window_shift_default
integration_windows_end = integration_windows_start + window_width_default

# then we define 2 possible edge cases
# the predicted integration window falls before the readout window
integration_before_readout = integration_windows_start < 0
# or after
integration_after_readout = integration_windows_end > (waveforms_to_repass.shape[1] - 1)

# If the resulting 5-samples window falls before the readout
# window we take the first 5 samples
window_width_before = 5
window_shift_before = 0

# in the case where the window is after, shift backward
window_width_after = 5
window_shift_after = 5

# and put them together:
window_widths = np.full(non_core_waveforms.shape[0], window_width_default)
window_widths[peak_before_readout] = window_width_before
window_widths[peak_after_readout] = window_width_after
window_shifts = np.full(non_core_waveforms.shape[0], window_shift_default)
window_shifts[peak_before_readout] = window_shift_before
window_shifts[peak_after_readout] = window_shift_after

# Now we can also (re)define the pathological predicted times because
# (we needed them to define the corrispective widths and shifts) set
# sample to 0 (beginning of the waveform) if predicted time falls before

predicted_peaks[predicted_peaks < 0] = 0
# If the resulting 5-samples window falls after the readout
# window we take the last 5 samples
window_width_after = 6
window_shift_after = 4

# put all values of widths and shifts for 2nd pass pixels together
window_widths = np.full(waveforms_to_repass.shape[0], window_width_default)
window_widths[integration_before_readout] = window_width_before
window_widths[integration_after_readout] = window_width_after
window_shifts = np.full(waveforms_to_repass.shape[0], window_shift_default)
window_shifts[integration_before_readout] = window_shift_before
window_shifts[integration_after_readout] = window_shift_after

# Now we have to (re)define the pathological predicted times for which
# - either the peak itself falls outside of the readout window
# - or is within the first or last 2 samples (so that at least 1 sample
# of the integration window is outside of the readout window)
# We place them at the first or last sample, so the special window
# widhts and shifts that we defined earlier put the integration window
# for these 2 cases either in the first 5 samples or the last

# set sample to 0 (beginning of the waveform)
# if predicted time falls before
# but also if it's so near the edge that the integration window falls
# outside
predicted_peaks[predicted_peaks < 2] = 0

# set sample to max-1 (first sample has index 0)
# if predicted time falls after
predicted_peaks[predicted_peaks > (waveforms.shape[1] - 1)] = (
waveforms.shape[1] - 1
predicted_peaks[predicted_peaks > (waveforms_to_repass.shape[1] - 3)] = (
waveforms_to_repass.shape[1] - 1
)

# re-calibrate non-core pixels using the fixed 5-samples window
charge_no_core, pulse_times_no_core = extract_around_peak(
non_core_waveforms,
# re-calibrate 2nd pass pixels using the fixed 5-samples window
reintegrated_charge, reestimated_pulse_times = extract_around_peak(
waveforms_to_repass,
predicted_peaks,
window_widths,
window_shifts,
Expand All @@ -1129,40 +1133,37 @@ def _apply_second_pass(
# now we compute 3 corrections for the default, before, and after cases:
correction = self._calculate_correction(
telid, window_width_default, window_shift_default
)[selected_gain_channel][non_core_pixels_mask]
)[selected_gain_channel][mask_2nd_pass]

correction_before = self._calculate_correction(
telid, window_width_before, window_shift_before
)[selected_gain_channel][non_core_pixels_mask]
)[selected_gain_channel][mask_2nd_pass]

correction_after = self._calculate_correction(
telid, window_width_after, window_shift_after
)[selected_gain_channel][non_core_pixels_mask]
)[selected_gain_channel][mask_2nd_pass]

correction[peak_before_readout] = correction_before[peak_before_readout]
correction[peak_after_readout] = correction_after[peak_after_readout]
correction[integration_before_readout] = correction_before[integration_before_readout]
correction[integration_after_readout] = correction_after[integration_after_readout]

charge_no_core *= correction
reintegrated_charge *= correction

# STEP 7

# Combine core and non-core pixels in the final output
# Combine in the final output with,
# - core pixels from the main cluster
# - rest of the pixels which have been passed a second time

# this is the biggest cluster from the cleaned image
# it contains the core pixels (which we leave untouched)
# plus possibly some non-core pixels
# Start from a copy of the 1st pass charge
charge_2ndpass = charge_1stpass.copy()
# Now we overwrite the charges of all non-core pixels in the camera
# plus all those pixels which didn't survive the preliminary
# cleaning.
# We apply also their corrections.
charge_2ndpass[non_core_pixels_mask] = charge_no_core
# Overwrite the charges of pixels marked for second pass
# leaving untouched the core pixels of the main island
# from the preliminary (cleaned) image
charge_2ndpass[mask_2nd_pass] = reintegrated_charge

# Same approach for the pulse times
pulse_time_2ndpass = pulse_time_1stpass # core + non-core pixels
pulse_time_2ndpass[
non_core_pixels_mask
] = pulse_times_no_core # non-core pixels
pulse_time_2ndpass = pulse_time_1stpass.copy()
pulse_time_2ndpass[mask_2nd_pass] = reestimated_pulse_times

return charge_2ndpass, pulse_time_2ndpass

Expand Down
Loading

0 comments on commit b4d1e29

Please sign in to comment.