sf-daq/dap
0
0
Fork 1
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

trailing whitespace

Browse Source
This commit is contained in:
augustin_s
2024-03-22 17:00:05 +01:00
parent c31e10a79f
commit 44af41b2f2
7 changed files with 120 additions and 120 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -11,7 +11,7 @@ clean:
rm -rf dist peakfinder8_extension.egg-info rm -rf dist peakfinder8_extension.egg-info
rm -rf peakfinder8/cython/peakfinder8_extension.cpp rm -rf peakfinder8/cython/peakfinder8_extension.cpp
peakfinder8/peakfinders.o: peakfinder8/cheetahmodules.h peakfinder8/peakfinders.cpp peakfinder8/peakfinders.h peakfinder8/peakfinders.o: peakfinder8/cheetahmodules.h peakfinder8/peakfinders.cpp peakfinder8/peakfinders.h
gcc -fPIC -I ./peakfinder8 -c peakfinder8/peakfinders.cpp -o peakfinder8/peakfinders.o gcc -fPIC -I ./peakfinder8 -c peakfinder8/peakfinders.cpp -o peakfinder8/peakfinders.o
peakfinder8/libpeakfinder8.so: peakfinder8/peakfinders.o peakfinder8/libpeakfinder8.so: peakfinder8/peakfinders.o

44
README.md
View File

@ -28,7 +28,7 @@ Clone and install dap
git clone https://gitlab.psi.ch/sf-daq/dap.git git clone https://gitlab.psi.ch/sf-daq/dap.git
cd dap cd dap
make install make install
``` ```
# Architecture # Architecture
@ -37,7 +37,7 @@ The dap architecture is designed for horizontal scalability, processing various
- Metadata-only results to the accumulator for storage. - Metadata-only results to the accumulator for storage.
## Worker ## Worker
Each worker runs independently and processes frames received via ZeroMQ. Before applying algorithms, it converts raw (ADC) detector values to energy using [jungfrau_utils](https://github.com/paulscherrerinstitute/jungfrau_utils). Each worker runs independently and processes frames received via ZeroMQ. Before applying algorithms, it converts raw (ADC) detector values to energy using [jungfrau_utils](https://github.com/paulscherrerinstitute/jungfrau_utils).
Input parameters: Input parameters:
```bash ```bash
python dap/worker.py --help python dap/worker.py --help
@ -65,7 +65,7 @@ options:
send to streamvis each of skip_frames_rate frames send to streamvis each of skip_frames_rate frames
``` ```
The number of required workers varies based on detector size and algorithm complexity. The number of required workers varies based on detector size and algorithm complexity.
Workers can be pinned to specific processor cores and distributed across multiple nodes. Workers can be pinned to specific processor cores and distributed across multiple nodes.
## Accumulator ## Accumulator
@ -87,9 +87,9 @@ options:
``` ```
# Implemented algorithms # Implemented algorithms
* **peakfinder Algorithm** * **peakfinder Algorithm**
This algorithm is based on peakfinder8 from the [cheetah package](https://www.desy.de/~barty/cheetah/Cheetah/Welcome.html). It identifies peaks as connected pixels exhibiting intensity above the background. The background is determined iteratively by radial averaging, excluding signal pixels. This algorithm is based on peakfinder8 from the [cheetah package](https://www.desy.de/~barty/cheetah/Cheetah/Welcome.html). It identifies peaks as connected pixels exhibiting intensity above the background. The background is determined iteratively by radial averaging, excluding signal pixels.
Input parameters: Input parameters:
@ -106,16 +106,16 @@ options:
* `'is_hit_frame': True/False` - Marks whether a frame qualifies as a hit based on the number of identified peaks exceeding the defined threshold. * `'is_hit_frame': True/False` - Marks whether a frame qualifies as a hit based on the number of identified peaks exceeding the defined threshold.
* **Radial Profile Integration** * **Radial Profile Integration**
This algorithm integrates pixel intensities radially based on defined parameters. This algorithm integrates pixel intensities radially based on defined parameters.
Input parameters: Input parameters:
* `'do_radial_integration': 1/0` - Indicates whether radial integration should occur within dap. * `'do_radial_integration': 1/0` - Indicates whether radial integration should occur within dap.
* `'beam_center_x/beam_center_y': float/float` - Specifies the beam center coordinates in the detector space. * `'beam_center_x/beam_center_y': float/float` - Specifies the beam center coordinates in the detector space.
* `'apply_threshold': 1/0` - Determines whether to apply a threshold to pixel intensities before radial integration. * `'apply_threshold': 1/0` - Determines whether to apply a threshold to pixel intensities before radial integration.
* `'radial_integration_silent_min/radial_integration_silent_max': float/float` - If both values are present, normalizes the radial integrated profile within this specified range. This is crucial for frame combination to eliminate variations in beam intensity across frames. * `'radial_integration_silent_min/radial_integration_silent_max': float/float` - If both values are present, normalizes the radial integrated profile within this specified range. This is crucial for frame combination to eliminate variations in beam intensity across frames.
Output of algorithm: Output of algorithm:
* `'radint_I': list[float]` - Represents the radial integrated profile in pixel coordinates. * `'radint_I': list[float]` - Represents the radial integrated profile in pixel coordinates.
* `'radint_q' : [float, float]` - Represents the minimum and maximum x-coordinate values considered during integration in pixel coordinates. * `'radint_q' : [float, float]` - Represents the minimum and maximum x-coordinate values considered during integration in pixel coordinates.
@ -139,7 +139,7 @@ options:
* `'roi_intensities': list[float]` - Sum of pixel intensities within the ROI. * `'roi_intensities': list[float]` - Sum of pixel intensities within the ROI.
* `'roi_intensities_normalised': list[float]` - Intensity within the defined ROI normalized to the ROI size (the count of active pixels within the ROI). * `'roi_intensities_normalised': list[float]` - Intensity within the defined ROI normalized to the ROI size (the count of active pixels within the ROI).
* `'roi_intensities_x': list(float, float)` - x1/x2 (left/right) x-coordinates of the ROI. * `'roi_intensities_x': list(float, float)` - x1/x2 (left/right) x-coordinates of the ROI.
* `'roi_intensities_proj_x': list(value)` - Projection onto the x-coordinate of pixel intensities (sum). * `'roi_intensities_proj_x': list(value)` - Projection onto the x-coordinate of pixel intensities (sum).
* **Detecting Frames with High intensity in Specific Regions** * **Detecting Frames with High intensity in Specific Regions**
@ -156,10 +156,10 @@ options:
* 0: if intensity in both ROIs falls below the respective thresholds * 0: if intensity in both ROIs falls below the respective thresholds
* 25: ROI1 has high energy but not ROI2 * 25: ROI1 has high energy but not ROI2
* 50: ROI2 has high energy but not ROI1 * 50: ROI2 has high energy but not ROI1
* 75: intensities in both ROIs exceed the thresholds * 75: intensities in both ROIs exceed the thresholds
* **Frame aggregation** * **Frame aggregation**
When dealing with faint signals that are challenging to discern in individual frames, dap offers the option to aggregate frames, combining them at the dap level before sending the resulting aggregate frame to visualization. It's important to note that this aggregation occurs independently for each worker. Thus, it's crucial to maintain a reasonable balance between the number of frames to aggregate and the number of active dap workers for a given detector. This aggregation process does not impact other algorithms, as they operate on individual frames (example: the threshold algorithm runs before frame aggregation). When dealing with faint signals that are challenging to discern in individual frames, dap offers the option to aggregate frames, combining them at the dap level before sending the resulting aggregate frame to visualization. It's important to note that this aggregation occurs independently for each worker. Thus, it's crucial to maintain a reasonable balance between the number of frames to aggregate and the number of active dap workers for a given detector. This aggregation process does not impact other algorithms, as they operate on individual frames (example: the threshold algorithm runs before frame aggregation).
Input parameters: Input parameters:
@ -167,22 +167,22 @@ options:
* `'aggregation_max': int` - Specifies the maximum number of frames to aggregate before transmitting to visualization. This value pertains to each worker. * `'aggregation_max': int` - Specifies the maximum number of frames to aggregate before transmitting to visualization. This value pertains to each worker.
Algorithm Output: Algorithm Output:
* `'aggregated_images': int` - Indicates the count of aggregated images. * `'aggregated_images': int` - Indicates the count of aggregated images.
* **Frame Tagging** * **Frame Tagging**
When event propagation is integrated into the detector, dap allows frames to be tagged accordingly, facilitating their categorization during visualization. When event propagation is integrated into the detector, dap allows frames to be tagged accordingly, facilitating their categorization during visualization.
Presently supported markers: Presently supported markers:
* `'laser_on': True/False` - Marks frames as "laser activated" when the darkshot event code is False and the laser event code is True. Otherwise, frames are labeled as "not laser activated". * `'laser_on': True/False` - Marks frames as "laser activated" when the darkshot event code is False and the laser event code is True. Otherwise, frames are labeled as "not laser activated".
* **Detection of Saturated Pixels** * **Detection of Saturated Pixels**
For every frame received by dap, an analysis is performed to ascertain the quantity and positions of saturated pixels. For every frame received by dap, an analysis is performed to ascertain the quantity and positions of saturated pixels.
Algorithm Output: Algorithm Output:
* `'saturated_pixels': int` - Number of saturated pixels within the frame. * `'saturated_pixels': int` - Number of saturated pixels within the frame.
* `'saturated_pixels_x/saturated_pixels_y': list[float]/list[float]` - Coordinates of the saturated pixels. * `'saturated_pixels_x/saturated_pixels_y': list[float]/list[float]` - Coordinates of the saturated pixels.
* **Transmitted Parameters from dap Input to Visualization** * **Transmitted Parameters from dap Input to Visualization**
@ -194,12 +194,12 @@ options:
* **Implement Additional Masking** * **Implement Additional Masking**
Sensitive algorithms, such as the peakfinder, often necessitate the exclusion of specific detector regions (like defective or module edge pixels). Presently, defining these detector segments requires manual hardcoding within the worker.py code. Sensitive algorithms, such as the peakfinder, often necessitate the exclusion of specific detector regions (like defective or module edge pixels). Presently, defining these detector segments requires manual hardcoding within the worker.py code.
Use the `'apply_additional_mask': 0/1` - Input flag to enable this functionality. Use the `'apply_additional_mask': 0/1` - Input flag to enable this functionality.
* **Filter based on pulse picker information** * **Filter based on pulse picker information**
If the event propagation capability is accessible for the detector and the pulse picker information is correctly configured for propagation, the filtration based on pulse picker information becomes feasible by using the If the event propagation capability is accessible for the detector and the pulse picker information is correctly configured for propagation, the filtration based on pulse picker information becomes feasible by using the
`'select_only_ppicker_events': 0/1` - Input flag. `'select_only_ppicker_events': 0/1` - Input flag.
## Input parameters (File) ## Input parameters (File)

6
dap/accumulator.py
View File

@ -14,14 +14,14 @@ def main():
parser.add_argument("--accumulator_port", default=13002, type=int, help="accumulator port") parser.add_argument("--accumulator_port", default=13002, type=int, help="accumulator port")
args = parser.parse_args() args = parser.parse_args()
FA_HOST_ACCUMULATE = args.accumulator FA_HOST_ACCUMULATE = args.accumulator
FA_PORT_ACCUMULATE = args.accumulator_port FA_PORT_ACCUMULATE = args.accumulator_port
zmq_context = zmq.Context(io_threads=4) zmq_context = zmq.Context(io_threads=4)
poller = zmq.Poller() poller = zmq.Poller()
# Accumulator # Accumulator
if True: if True:
accumulator_socket = zmq_context.socket(zmq.PULL) accumulator_socket = zmq_context.socket(zmq.PULL)
accumulator_socket.bind('tcp://*:%s' % FA_PORT_ACCUMULATE ) accumulator_socket.bind('tcp://*:%s' % FA_PORT_ACCUMULATE )
@ -42,7 +42,7 @@ def main():
if accumulator_socket in events: if accumulator_socket in events:
results = accumulator_socket.recv_json(flags) results = accumulator_socket.recv_json(flags)
n_frames_received += 1 n_frames_received += 1
pulse_id = results.get('pulse_id', 0) pulse_id = results.get('pulse_id', 0)
run_name = str(pulse_id//10000*10000) run_name = str(pulse_id//10000*10000)
detector = results.get('detector_name', "") detector = results.get('detector_name', "")

44
dap/worker.py
View File

@ -23,7 +23,7 @@ def radial_profile(data, r, nr, keep_pixels=None):
else: else:
tbin = np.bincount(r, data.ravel()) tbin = np.bincount(r, data.ravel())
radialprofile = tbin / nr radialprofile = tbin / nr
return radialprofile return radialprofile
def prepare_radial_profile(data, center, keep_pixels=None): def prepare_radial_profile(data, center, keep_pixels=None):
y, x = np.indices((data.shape)) y, x = np.indices((data.shape))
@ -54,7 +54,7 @@ def main():
parser.add_argument("--skip_frames_rate", default=1, type=int, help="send to streamvis each of skip_frames_rate frames") parser.add_argument("--skip_frames_rate", default=1, type=int, help="send to streamvis each of skip_frames_rate frames")
args = parser.parse_args() args = parser.parse_args()
if args.backend: if args.backend:
BACKEND_ADDRESS = args.backend BACKEND_ADDRESS = args.backend
else: else:
@ -87,7 +87,7 @@ def main():
worker = 1 worker = 1
# receive from backend: # receive from backend:
backend_socket = zmq_context.socket(zmq.PULL) backend_socket = zmq_context.socket(zmq.PULL)
backend_socket.connect(BACKEND_ADDRESS) backend_socket.connect(BACKEND_ADDRESS)
poller.register(backend_socket, zmq.POLLIN) poller.register(backend_socket, zmq.POLLIN)
@ -140,7 +140,7 @@ def main():
metadata = backend_socket.recv_json(flags) metadata = backend_socket.recv_json(flags)
image = backend_socket.recv(flags, copy=False, track=False) image = backend_socket.recv(flags, copy=False, track=False)
image = np.frombuffer(image, dtype=metadata['type']).reshape(metadata['shape']) image = np.frombuffer(image, dtype=metadata['type']).reshape(metadata['shape'])
results = copy(metadata) results = copy(metadata)
if results['shape'][0] == 2 and results['shape'][1] == 2: if results['shape'][0] == 2 and results['shape'][1] == 2:
@ -160,7 +160,7 @@ def main():
event_darkshot = bool((daq_rec>>17)&1) event_darkshot = bool((daq_rec>>17)&1)
event_fel = bool((daq_rec>>18)&1) event_fel = bool((daq_rec>>18)&1)
event_ppicker = bool((daq_rec>>19)&1) event_ppicker = bool((daq_rec>>19)&1)
# #
if not event_darkshot: if not event_darkshot:
results['laser_on'] = event_laser results['laser_on'] = event_laser
@ -174,7 +174,7 @@ def main():
# if event_ppicker: # if event_ppicker:
# results['number_of_spots'] = 50 # results['number_of_spots'] = 50
# results['is_hit_frame'] = True # results['is_hit_frame'] = True
double_pixels = results.get('double_pixels', "mask") double_pixels = results.get('double_pixels', "mask")
pedestal_file_name = metadata.get("pedestal_name", None) pedestal_file_name = metadata.get("pedestal_name", None)
@ -190,12 +190,12 @@ def main():
# pedestal file is not in stream, skip this frame # pedestal file is not in stream, skip this frame
if ju_stream_adapter.handler.pedestal_file is None or ju_stream_adapter.handler.pedestal_file == "": if ju_stream_adapter.handler.pedestal_file is None or ju_stream_adapter.handler.pedestal_file == "":
continue continue
data=np.ascontiguousarray(data) data=np.ascontiguousarray(data)
# starting from ju 3.3.1 pedestal file is cached in library, re-calculated only if parameters(and/or pedestal file) are changed # starting from ju 3.3.1 pedestal file is cached in library, re-calculated only if parameters(and/or pedestal file) are changed
id_pixel_mask_1 = id(pixel_mask_corrected) id_pixel_mask_1 = id(pixel_mask_corrected)
pixel_mask_corrected=ju_stream_adapter.handler.get_pixel_mask(geometry=True, gap_pixels=True, double_pixels=double_pixels) pixel_mask_corrected=ju_stream_adapter.handler.get_pixel_mask(geometry=True, gap_pixels=True, double_pixels=double_pixels)
id_pixel_mask_2 = id(pixel_mask_corrected) id_pixel_mask_2 = id(pixel_mask_corrected)
if id_pixel_mask_1 != id_pixel_mask_2: if id_pixel_mask_1 != id_pixel_mask_2:
@ -206,14 +206,14 @@ def main():
pixel_mask_pf = np.ascontiguousarray(pixel_mask_corrected).astype(np.int8, copy=False) pixel_mask_pf = np.ascontiguousarray(pixel_mask_corrected).astype(np.int8, copy=False)
else: else:
pixel_mask_pf = None pixel_mask_pf = None
# #
disabled_modules = results.get("disabled_modules", []) disabled_modules = results.get("disabled_modules", [])
# add additional mask at the edge of modules for JF06T08 # add additional mask at the edge of modules for JF06T08
apply_additional_mask = (results.get("apply_additional_mask", 0) == 1) apply_additional_mask = (results.get("apply_additional_mask", 0) == 1)
if detector == "JF06T08V04" and apply_additional_mask: if detector == "JF06T08V04" and apply_additional_mask:
# edge pixels # edge pixels
pixel_mask_pf[67:1097,1063] = 0 pixel_mask_pf[67:1097,1063] = 0
pixel_mask_pf[0:1030, 1100] = 0 pixel_mask_pf[0:1030, 1100] = 0
@ -232,7 +232,7 @@ def main():
pixel_mask_pf[67:1097, 550] = 0 pixel_mask_pf[67:1097, 550] = 0
pixel_mask_pf[0:1030, 1650] = 0 pixel_mask_pf[0:1030, 1650] = 0
pixel_mask_pf[0:1030, 1613] = 0 pixel_mask_pf[0:1030, 1613] = 0
pixel_mask_pf[1106, 68:582] = 0 pixel_mask_pf[1106, 68:582] = 0
@ -267,7 +267,7 @@ def main():
if detector == "JF17T16V01" and apply_additional_mask: if detector == "JF17T16V01" and apply_additional_mask:
# mask module 11 # mask module 11
pixel_mask_pf[2619:3333,1577:2607] = 0 pixel_mask_pf[2619:3333,1577:2607] = 0
if pixel_mask_corrected is not None: if pixel_mask_corrected is not None:
data_s = copy(image) data_s = copy(image)
saturated_pixels_coordinates = ju_stream_adapter.handler.get_saturated_pixels(data_s, mask=True, geometry=True, gap_pixels=True, double_pixels=double_pixels) saturated_pixels_coordinates = ju_stream_adapter.handler.get_saturated_pixels(data_s, mask=True, geometry=True, gap_pixels=True, double_pixels=double_pixels)
@ -277,7 +277,7 @@ def main():
# pump probe analysis # pump probe analysis
do_radial_integration = results.get("do_radial_integration", 0) do_radial_integration = results.get("do_radial_integration", 0)
if (do_radial_integration != 0): if (do_radial_integration != 0):
data_copy_1 = np.copy(data) data_copy_1 = np.copy(data)
@ -290,9 +290,9 @@ def main():
if r_radial_integration is None: if r_radial_integration is None:
r_radial_integration, nr_radial_integration = prepare_radial_profile(data_copy_1, center_radial_integration, keep_pixels) r_radial_integration, nr_radial_integration = prepare_radial_profile(data_copy_1, center_radial_integration, keep_pixels)
r_min_max = [int(np.min(r_radial_integration)), int(np.max(r_radial_integration))+1] r_min_max = [int(np.min(r_radial_integration)), int(np.max(r_radial_integration))+1]
apply_threshold = results.get('apply_threshold', 0)
apply_threshold = results.get('apply_threshold', 0)
if (apply_threshold != 0) and all ( k in results for k in ('threshold_min', 'threshold_max')): if (apply_threshold != 0) and all ( k in results for k in ('threshold_min', 'threshold_max')):
threshold_min = float(results['threshold_min']) threshold_min = float(results['threshold_min'])
@ -305,11 +305,11 @@ def main():
silent_region_min = results.get("radial_integration_silent_min", None) silent_region_min = results.get("radial_integration_silent_min", None)
silent_region_max = results.get("radial_integration_silent_max", None) silent_region_max = results.get("radial_integration_silent_max", None)
if ( silent_region_min is not None and silent_region_max is not None and if ( silent_region_min is not None and silent_region_max is not None and
silent_region_max > silent_region_min and silent_region_max > silent_region_min and
silent_region_min > r_min_max[0] and silent_region_max < r_min_max[1] ): silent_region_min > r_min_max[0] and silent_region_max < r_min_max[1] ):
integral_silent_region = np.sum(rp[silent_region_min:silent_region_max]) integral_silent_region = np.sum(rp[silent_region_min:silent_region_max])
rp = rp/integral_silent_region rp = rp/integral_silent_region
results['radint_normalised'] = [silent_region_min, silent_region_max] results['radint_normalised'] = [silent_region_min, silent_region_max]
@ -352,7 +352,7 @@ def main():
for iRoi in range(len(roi_x1)): for iRoi in range(len(roi_x1)):
data_roi = np.copy(d[roi_y1[iRoi]:roi_y2[iRoi],roi_x1[iRoi]:roi_x2[iRoi]]) data_roi = np.copy(d[roi_y1[iRoi]:roi_y2[iRoi],roi_x1[iRoi]:roi_x2[iRoi]])
roi_results[iRoi] = np.nansum(data_roi) roi_results[iRoi] = np.nansum(data_roi)
if threshold_value_choice == "NaN": if threshold_value_choice == "NaN":
roi_results_normalised[iRoi] = roi_results[iRoi]/((roi_y2[iRoi]-roi_y1[iRoi])*(roi_x2[iRoi]-roi_x1[iRoi])) roi_results_normalised[iRoi] = roi_results[iRoi]/((roi_y2[iRoi]-roi_y1[iRoi])*(roi_x2[iRoi]-roi_x1[iRoi]))
else: else:
@ -377,7 +377,7 @@ def main():
if results['roi_intensities_normalised'][1] >= results['spi_limit'][1]: if results['roi_intensities_normalised'][1] >= results['spi_limit'][1]:
number_of_spots += 50 number_of_spots += 50
results['number_of_spots'] = number_of_spots results['number_of_spots'] = number_of_spots
if number_of_spots > 0: if number_of_spots > 0:
results['is_hit_frame'] = True results['is_hit_frame'] = True

2
peakfinder8/cheetahmodules.h
View File

@ -16,7 +16,7 @@ int killNearbyPeaks(tPeakList*, float );
#include <stdint.h> #include <stdint.h>
// from detectorObject.h file // from detectorObject.h file
//-------------------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------------------
/* /*
* Bits for pixel masks * Bits for pixel masks
* Oriented along conventions defined for CXI file format ( https://github.com/FilipeMaia/CXI/raw/master/cxi_file_format.pdf ) * Oriented along conventions defined for CXI file format ( https://github.com/FilipeMaia/CXI/raw/master/cxi_file_format.pdf )

140
peakfinder8/peakfinders.cpp
View File

@ -21,14 +21,14 @@
*/ */
void allocatePeakList(tPeakList *peak, long NpeaksMax) { void allocatePeakList(tPeakList *peak, long NpeaksMax) {
peak->nPeaks = 0; peak->nPeaks = 0;
peak->nPeaks_max = NpeaksMax; peak->nPeaks_max = NpeaksMax;
peak->nHot = 0; peak->nHot = 0;
peak->peakResolution = 0; peak->peakResolution = 0;
peak->peakResolutionA = 0; peak->peakResolutionA = 0;
peak->peakDensity = 0; peak->peakDensity = 0;
peak->peakNpix = 0; peak->peakNpix = 0;
peak->peakTotal = 0; peak->peakTotal = 0;
peak->peak_maxintensity = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_maxintensity = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_totalintensity = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_totalintensity = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_sigma = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_sigma = (float *) calloc(NpeaksMax, sizeof(float));
@ -41,7 +41,7 @@ void allocatePeakList(tPeakList *peak, long NpeaksMax) {
peak->peak_com_y_assembled = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_com_y_assembled = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_com_r_assembled = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_com_r_assembled = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_com_q = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_com_q = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_com_res = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_com_res = (float *) calloc(NpeaksMax, sizeof(float));
peak->memoryAllocated = 1; peak->memoryAllocated = 1;
} }
@ -75,20 +75,20 @@ void freePeakList(tPeakList peak) {
* Anton Barty * Anton Barty
*/ */
int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long asic_nx, long asic_ny, long nasics_x, long nasics_y, float ADCthresh, float hitfinderMinSNR, long hitfinderMinPixCount, long hitfinderMaxPixCount, long hitfinderLocalBGRadius) { int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long asic_nx, long asic_ny, long nasics_x, long nasics_y, float ADCthresh, float hitfinderMinSNR, long hitfinderMinPixCount, long hitfinderMaxPixCount, long hitfinderLocalBGRadius) {
// Derived values // Derived values
long pix_nx = asic_nx*nasics_x; long pix_nx = asic_nx*nasics_x;
long pix_ny = asic_ny*nasics_y; long pix_ny = asic_ny*nasics_y;
long pix_nn = pix_nx*pix_ny; long pix_nn = pix_nx*pix_ny;
//long asic_nn = asic_nx*asic_ny; //long asic_nn = asic_nx*asic_ny;
long hitfinderNpeaksMax = peaklist->nPeaks_max; long hitfinderNpeaksMax = peaklist->nPeaks_max;
peaklist->nPeaks = 0; peaklist->nPeaks = 0;
peaklist->peakNpix = 0; peaklist->peakNpix = 0;
peaklist->peakTotal = 0; peaklist->peakTotal = 0;
// Variables for this hitfinder // Variables for this hitfinder
long nat = 0; long nat = 0;
long lastnat = 0; long lastnat = 0;
@ -111,28 +111,28 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
long fs, ss; long fs, ss;
float com_x, com_y, com_e; float com_x, com_y, com_e;
float thisADCthresh; float thisADCthresh;
nat = 0; nat = 0;
//counter = 0; //counter = 0;
total = 0.0; total = 0.0;
snr=0; snr=0;
maxI = 0; maxI = 0;
/* /*
* Create a buffer for image data so we don't nuke the main image by mistake * Create a buffer for image data so we don't nuke the main image by mistake
*/ */
float *temp = (float*) malloc(pix_nn*sizeof(float)); float *temp = (float*) malloc(pix_nn*sizeof(float));
memcpy(temp, data, pix_nn*sizeof(float)); memcpy(temp, data, pix_nn*sizeof(float));
/* /*
* Apply mask (multiply data by 0 to ignore regions - this makes data below threshold for peak finding) * Apply mask (multiply data by 0 to ignore regions - this makes data below threshold for peak finding)
*/ */
for(long i=0;i<pix_nn;i++){ for(long i=0;i<pix_nn;i++){
temp[i] *= mask[i]; temp[i] *= mask[i];
} }
/* /*
* Determine noise and offset as a funciton of radius * Determine noise and offset as a funciton of radius
*/ */
@ -140,7 +140,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
long lminr, lmaxr; long lminr, lmaxr;
fminr = 1e9; fminr = 1e9;
fmaxr = -1e9; fmaxr = -1e9;
// Figure out radius bounds // Figure out radius bounds
for(long i=0;i<pix_nn;i++){ for(long i=0;i<pix_nn;i++){
if (pix_r[i] > fmaxr) if (pix_r[i] > fmaxr)
@ -150,16 +150,16 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
} }
lmaxr = (int)ceil(fmaxr)+1; lmaxr = (int)ceil(fmaxr)+1;
lminr = 0; lminr = 0;
// Allocate and zero arrays // Allocate and zero arrays
float *rsigma = (float*) malloc(lmaxr*sizeof(float)); float *rsigma = (float*) malloc(lmaxr*sizeof(float));
float *roffset = (float*) malloc(lmaxr*sizeof(float)); float *roffset = (float*) malloc(lmaxr*sizeof(float));
long *rcount = (long*) malloc(lmaxr*sizeof(long)); long *rcount = (long*) malloc(lmaxr*sizeof(long));
float *rthreshold = (float*) malloc(lmaxr*sizeof(float)); float *rthreshold = (float*) malloc(lmaxr*sizeof(float));
long *peakpixels = (long *) calloc(hitfinderMaxPixCount, sizeof(long)); long *peakpixels = (long *) calloc(hitfinderMaxPixCount, sizeof(long));
char *peakpixel = (char *) calloc(pix_nn, sizeof(char)); char *peakpixel = (char *) calloc(pix_nn, sizeof(char));
char *rthreshold_changed = (char *) malloc(lmaxr*sizeof(char)); char *rthreshold_changed = (char *) malloc(lmaxr*sizeof(char));
int *pix_rint = (int *) malloc(pix_nn*sizeof(int)); int *pix_rint = (int *) malloc(pix_nn*sizeof(int));
@ -177,7 +177,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
pixels_check[i] = i; pixels_check[i] = i;
} }
long n_pixels_check = pix_nn; long n_pixels_check = pix_nn;
// Compute sigma and average of data values at each radius // Compute sigma and average of data values at each radius
// From this, compute the ADC threshold to be applied at each radius // From this, compute the ADC threshold to be applied at each radius
// Iterate a few times to reduce the effect of positive outliers (ie: peaks) // Iterate a few times to reduce the effect of positive outliers (ie: peaks)
@ -215,7 +215,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
rcount[thisr] += 1; rcount[thisr] += 1;
} }
pixels_check[new_pixels_check] = i; pixels_check[new_pixels_check] = i;
new_pixels_check++; new_pixels_check++;
} }
} }
} }
@ -249,7 +249,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
} }
} }
} }
com_x=0; com_x=0;
com_y=0; com_y=0;
@ -257,24 +257,24 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
for(long mj=0; mj<nasics_y; mj++){ for(long mj=0; mj<nasics_y; mj++){
for(long mi=0; mi<nasics_x; mi++){ for(long mi=0; mi<nasics_x; mi++){
// Loop over pixels within a module // Loop over pixels within a module
for(long j=1; j<asic_ny-1; j++){ for(long j=1; j<asic_ny-1; j++){
for(long i=1; i<asic_nx-1; i++){ for(long i=1; i<asic_nx-1; i++){
ss = (j+mj*asic_ny)*pix_nx; ss = (j+mj*asic_ny)*pix_nx;
fs = i+mi*asic_nx; fs = i+mi*asic_nx;
e = ss + fs; e = ss + fs;
if(e > pix_nn) { if(e > pix_nn) {
printf("Array bounds error: e=%li\n",e); printf("Array bounds error: e=%li\n",e);
exit(1); exit(1);
} }
thisr = pix_rint[e]; thisr = pix_rint[e];
thisADCthresh = rthreshold[thisr]; thisADCthresh = rthreshold[thisr];
if(temp[e] > thisADCthresh && peakpixel[e] == 0){ if(temp[e] > thisADCthresh && peakpixel[e] == 0){
// This might be the start of a new peak - start searching // This might be the start of a new peak - start searching
inx[0] = i; inx[0] = i;
@ -287,10 +287,10 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
maxIraw = 0; maxIraw = 0;
peak_com_x = 0; peak_com_x = 0;
peak_com_y = 0; peak_com_y = 0;
// Keep looping until the pixel count within this peak does not change // Keep looping until the pixel count within this peak does not change
do { do {
lastnat = nat; lastnat = nat;
// Loop through points known to be within this peak // Loop through points known to be within this peak
for(long p=0; p<nat; p++){ for(long p=0; p<nat; p++){
@ -305,20 +305,20 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
continue; continue;
if((iny[p]+search_y[k]) >= asic_ny) if((iny[p]+search_y[k]) >= asic_ny)
continue; continue;
// Neighbour point in big array // Neighbour point in big array
thisx = inx[p]+search_x[k]+mi*asic_nx; thisx = inx[p]+search_x[k]+mi*asic_nx;
thisy = iny[p]+search_y[k]+mj*asic_ny; thisy = iny[p]+search_y[k]+mj*asic_ny;
e = thisx + thisy*pix_nx; e = thisx + thisy*pix_nx;
//if(e < 0 || e >= pix_nn){ //if(e < 0 || e >= pix_nn){
// printf("Array bounds error: e=%i\n",e); // printf("Array bounds error: e=%i\n",e);
// continue; // continue;
//} //}
thisr = pix_rint[e]; thisr = pix_rint[e];
thisADCthresh = rthreshold[thisr]; thisADCthresh = rthreshold[thisr];
// Above threshold? // Above threshold?
if(temp[e] > thisADCthresh && peakpixel[e] == 0 && mask[e] != 0){ if(temp[e] > thisADCthresh && peakpixel[e] == 0 && mask[e] != 0){
//if(nat < 0 || nat >= global->pix_nn) { //if(nat < 0 || nat >= global->pix_nn) {
@ -340,20 +340,20 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
maxI = thisI; maxI = thisI;
if (thisI > maxIraw) if (thisI > maxIraw)
maxIraw = temp[e]; maxIraw = temp[e];
nat++; nat++;
} }
} }
} }
} while(lastnat != nat); } while(lastnat != nat);
// Too many or too few pixels means ignore this 'peak'; move on now // Too many or too few pixels means ignore this 'peak'; move on now
if(nat<hitfinderMinPixCount || nat>hitfinderMaxPixCount) { if(nat<hitfinderMinPixCount || nat>hitfinderMaxPixCount) {
continue; continue;
} }
/* /*
* Calculate center of mass for this peak from initial peak search * Calculate center of mass for this peak from initial peak search
*/ */
@ -363,8 +363,8 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
long com_xi = lrint(com_x) - mi*asic_nx; long com_xi = lrint(com_x) - mi*asic_nx;
long com_yi = lrint(com_y) - mj*asic_ny; long com_yi = lrint(com_y) - mj*asic_ny;
/* /*
* Calculate the local signal-to-noise ratio and local background in an annulus around this peak * Calculate the local signal-to-noise ratio and local background in an annulus around this peak
* (excluding pixels which look like they might be part of another peak) * (excluding pixels which look like they might be part of another peak)
@ -372,7 +372,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
float localSigma=0; float localSigma=0;
float localOffset=0; float localOffset=0;
long ringWidth = 2*hitfinderLocalBGRadius; long ringWidth = 2*hitfinderLocalBGRadius;
float sumI = 0; float sumI = 0;
float sumIsquared = 0; float sumIsquared = 0;
long np_sigma = 0; long np_sigma = 0;
@ -380,10 +380,10 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
float fbgr; float fbgr;
float backgroundMaxI=0; float backgroundMaxI=0;
float fBackgroundThresh=0; float fBackgroundThresh=0;
for(long bj=-ringWidth; bj<ringWidth; bj++){ for(long bj=-ringWidth; bj<ringWidth; bj++){
for(long bi=-ringWidth; bi<ringWidth; bi++){ for(long bi=-ringWidth; bi<ringWidth; bi++){
// Within-ASIC check // Within-ASIC check
if((com_xi+bi) < 0) if((com_xi+bi) < 0)
continue; continue;
@ -393,28 +393,28 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
continue; continue;
if((com_yi+bj) >= asic_ny) if((com_yi+bj) >= asic_ny)
continue; continue;
// Within outer ring check // Within outer ring check
fbgr = sqrt( bi*bi + bj*bj ); fbgr = sqrt( bi*bi + bj*bj );
if( fbgr > ringWidth )// || fbgr <= hitfinderLocalBGRadius ) // || fbgr > hitfinderLocalBGRadius) if( fbgr > ringWidth )// || fbgr <= hitfinderLocalBGRadius ) // || fbgr > hitfinderLocalBGRadius)
continue; continue;
// Position of this point in data stream // Position of this point in data stream
thisx = com_xi + bi + mi*asic_nx; thisx = com_xi + bi + mi*asic_nx;
thisy = com_yi + bj + mj*asic_ny; thisy = com_yi + bj + mj*asic_ny;
e = thisx + thisy*pix_nx; e = thisx + thisy*pix_nx;
thisr = pix_rint[e]; thisr = pix_rint[e];
thisADCthresh = rthreshold[thisr]; thisADCthresh = rthreshold[thisr];
// Intensity above background // Intensity above background
thisI = temp[e]; thisI = temp[e];
// If above ADC threshold, this could be part of another peak // If above ADC threshold, this could be part of another peak
//if (temp[e] > thisADCthresh) //if (temp[e] > thisADCthresh)
// continue; // continue;
// Keep track of value and value-squared for offset and sigma calculation // Keep track of value and value-squared for offset and sigma calculation
// if(peakpixel[e] == 0 && mask[e] != 0) { // if(peakpixel[e] == 0 && mask[e] != 0) {
if(temp[e] < thisADCthresh && peakpixel[e] == 0 && mask[e] != 0) { if(temp[e] < thisADCthresh && peakpixel[e] == 0 && mask[e] != 0) {
@ -428,7 +428,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
np_counted += 1; np_counted += 1;
} }
} }
// Calculate local background and standard deviation // Calculate local background and standard deviation
if (np_sigma != 0) { if (np_sigma != 0) {
localOffset = sumI/np_sigma; localOffset = sumI/np_sigma;
@ -438,8 +438,8 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
localOffset = roffset[pix_rint[lrint(com_e)]]; localOffset = roffset[pix_rint[lrint(com_e)]];
localSigma = 0.01; localSigma = 0.01;
} }
/* /*
* Re-integrate (and re-centroid) peak using local background estimates * Re-integrate (and re-centroid) peak using local background estimates
*/ */
@ -453,7 +453,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
e = peakpixels[counter]; e = peakpixels[counter];
thisIraw = temp[e]; thisIraw = temp[e];
thisI = thisIraw - localOffset; thisI = thisIraw - localOffset;
totI += thisI; totI += thisI;
totIraw += thisIraw; totIraw += thisIraw;
@ -473,7 +473,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
com_y = peak_com_y/fabs(totI); com_y = peak_com_y/fabs(totI);
com_e = lrint(com_x) + lrint(com_y)*pix_nx; com_e = lrint(com_x) + lrint(com_y)*pix_nx;
/* /*
* Calculate signal-to-noise and apply SNR criteria * Calculate signal-to-noise and apply SNR criteria
@ -482,11 +482,11 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
//snr = (float) (maxI)/localSigma; //snr = (float) (maxI)/localSigma;
//snr = (float) (totIraw-nat*localOffset)/localSigma; //snr = (float) (totIraw-nat*localOffset)/localSigma;
//snr = (float) (maxIraw-localOffset)/localSigma; //snr = (float) (maxIraw-localOffset)/localSigma;
// The more pixels there are in the peak, the more relaxed we are about this criterion // The more pixels there are in the peak, the more relaxed we are about this criterion
if( snr < hitfinderMinSNR ) // - nat +hitfinderMinPixCount if( snr < hitfinderMinSNR ) // - nat +hitfinderMinPixCount
continue; continue;
// Is the maximum intensity in the peak enough above intensity in background region to be a peak and not noise? // Is the maximum intensity in the peak enough above intensity in background region to be a peak and not noise?
// The more pixels there are in the peak, the more relaxed we are about this criterion // The more pixels there are in the peak, the more relaxed we are about this criterion
//fBackgroundThresh = hitfinderMinSNR - nat; //fBackgroundThresh = hitfinderMinSNR - nat;
@ -495,24 +495,24 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
fBackgroundThresh *= (backgroundMaxI-localOffset); fBackgroundThresh *= (backgroundMaxI-localOffset);
if( maxI < fBackgroundThresh) if( maxI < fBackgroundThresh)
continue; continue;
// This is a peak? If so, add info to peak list // This is a peak? If so, add info to peak list
if(nat>=hitfinderMinPixCount && nat<=hitfinderMaxPixCount ) { if(nat>=hitfinderMinPixCount && nat<=hitfinderMaxPixCount ) {
// This CAN happen! // This CAN happen!
if(totI == 0) if(totI == 0)
continue; continue;
//com_x = peak_com_x/fabs(totI); //com_x = peak_com_x/fabs(totI);
//com_y = peak_com_y/fabs(totI); //com_y = peak_com_y/fabs(totI);
e = lrint(com_x) + lrint(com_y)*pix_nx; e = lrint(com_x) + lrint(com_y)*pix_nx;
if(e < 0 || e >= pix_nn){ if(e < 0 || e >= pix_nn){
printf("Array bounds error: e=%ld\n",e); printf("Array bounds error: e=%ld\n",e);
continue; continue;
} }
// Remember peak information // Remember peak information
if (peakCounter < hitfinderNpeaksMax) { if (peakCounter < hitfinderNpeaksMax) {
peaklist->peakNpix += nat; peaklist->peakNpix += nat;
@ -539,28 +539,28 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
} }
//END: ; //END: ;
free(temp); free(temp);
free(inx); free(inx);
free(iny); free(iny);
free(peakpixel); free(peakpixel);
free(peakpixels); free(peakpixels);
free(roffset); free(roffset);
free(rsigma); free(rsigma);
free(rcount); free(rcount);
free(rthreshold); free(rthreshold);
free(pix_rint); free(pix_rint);
free(pixels_check); free(pixels_check);
free(rthreshold_changed); free(rthreshold_changed);
peaklist->nPeaks = peakCounter; peaklist->nPeaks = peakCounter;
return(peaklist->nPeaks); return(peaklist->nPeaks);
/*************************************************/ /*************************************************/
} }

2
peakfinder8/peakfinders.h
View File

@ -22,7 +22,7 @@ public:
float peakTotal; // Total integrated intensity in peaks float peakTotal; // Total integrated intensity in peaks
int memoryAllocated; int memoryAllocated;
long nPeaks_max; long nPeaks_max;
float *peak_maxintensity; // Maximum intensity in peak float *peak_maxintensity; // Maximum intensity in peak
float *peak_totalintensity; // Integrated intensity in peak float *peak_totalintensity; // Integrated intensity in peak
float *peak_sigma; // Signal-to-noise ratio of peak float *peak_sigma; // Signal-to-noise ratio of peak