From 942eea9808c0b8ffaa25e199903b4e8ea7e31799 Mon Sep 17 00:00:00 2001 From: Dmitry Ozerov Date: Thu, 28 Sep 2023 16:10:55 +0200 Subject: [PATCH] first commit --- Makefile | 23 + README.md | 99 +--- ap/accumulator.py | 66 +++ ap/worker.py | 501 ++++++++++++++++ peakfinder8/.gitignore | 5 + peakfinder8/cheetahmodules.h | 46 ++ peakfinder8/cython/peakfinder8.hh | 1 + peakfinder8/cython/peakfinder8_extension.pyx | 90 +++ peakfinder8/peakfinders.cpp | 566 +++++++++++++++++++ peakfinder8/peakfinders.h | 54 ++ setup.py | 15 + 11 files changed, 1388 insertions(+), 78 deletions(-) create mode 100644 Makefile create mode 100644 ap/accumulator.py create mode 100644 ap/worker.py create mode 100644 peakfinder8/.gitignore create mode 100644 peakfinder8/cheetahmodules.h create mode 100644 peakfinder8/cython/peakfinder8.hh create mode 100644 peakfinder8/cython/peakfinder8_extension.pyx create mode 100644 peakfinder8/peakfinders.cpp create mode 100644 peakfinder8/peakfinders.h create mode 100644 setup.py diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..e1e9aa2 --- /dev/null +++ b/Makefile @@ -0,0 +1,23 @@ + +all: peakfinder8_extension.so + +peakfinder8_extension.so: peakfinder8/libpeakfinder8.so + python setup.py build_ext --inplace + +clean: + rm -rf build + rm -rf peakfinder8/libpeakfinder8.so peakfinder8/peakfinders.o + rm -rf peakfinder8_extension.*.so + rm -rf dist peakfinder8_extension.egg-info + rm -rf peakfinder8/cython/peakfinder8_extension.cpp + +peakfinder8/peakfinders.o: peakfinder8/cheetahmodules.h peakfinder8/peakfinders.cpp peakfinder8/peakfinders.h + gcc -fPIC -I ./peakfinder8 -c peakfinder8/peakfinders.cpp -o peakfinder8/peakfinders.o + +peakfinder8/libpeakfinder8.so: peakfinder8/peakfinders.o + gcc -shared -o peakfinder8/libpeakfinder8.so peakfinder8/peakfinders.o -lc + +install: all + cp peakfinder8/libpeakfinder8.so ${CONDA_PREFIX}/lib/. + python setup.py install + diff --git a/README.md b/README.md index 72f32f8..f57a76e 100644 --- a/README.md +++ b/README.md @@ -3,91 +3,34 @@ Detector Analysis Pipeline Runs on detector data stream produced by sf-daq -## Getting started - -To make it easy for you to get started with GitLab, here's a list of recommended next steps. - -Already a pro? Just edit this README.md and make it your own. Want to make it easy? [Use the template at the bottom](#editing-this-readme)! - -## Add your files - -- [ ] [Create](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#create-a-file) or [upload](https://docs.gitlab.com/ee/user/project/repository/web_editor.html#upload-a-file) files -- [ ] [Add files using the command line](https://docs.gitlab.com/ee/gitlab-basics/add-file.html#add-a-file-using-the-command-line) or push an existing Git repository with the following command: - -``` -cd existing_repo -git remote add origin https://gitlab.psi.ch/sf-daq/dap.git -git branch -M main -git push -uf origin main -``` - -## Integrate with your tools - -- [ ] [Set up project integrations](https://gitlab.psi.ch/sf-daq/dap/-/settings/integrations) - -## Collaborate with your team - -- [ ] [Invite team members and collaborators](https://docs.gitlab.com/ee/user/project/members/) -- [ ] [Create a new merge request](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html) -- [ ] [Automatically close issues from merge requests](https://docs.gitlab.com/ee/user/project/issues/managing_issues.html#closing-issues-automatically) -- [ ] [Enable merge request approvals](https://docs.gitlab.com/ee/user/project/merge_requests/approvals/) -- [ ] [Automatically merge when pipeline succeeds](https://docs.gitlab.com/ee/user/project/merge_requests/merge_when_pipeline_succeeds.html) - -## Test and Deploy - -Use the built-in continuous integration in GitLab. - -- [ ] [Get started with GitLab CI/CD](https://docs.gitlab.com/ee/ci/quick_start/index.html) -- [ ] [Analyze your code for known vulnerabilities with Static Application Security Testing(SAST)](https://docs.gitlab.com/ee/user/application_security/sast/) -- [ ] [Deploy to Kubernetes, Amazon EC2, or Amazon ECS using Auto Deploy](https://docs.gitlab.com/ee/topics/autodevops/requirements.html) -- [ ] [Use pull-based deployments for improved Kubernetes management](https://docs.gitlab.com/ee/user/clusters/agent/) -- [ ] [Set up protected environments](https://docs.gitlab.com/ee/ci/environments/protected_environments.html) - -*** - -# Editing this README - -When you're ready to make this README your own, just edit this file and use the handy template below (or feel free to structure it however you want - this is just a starting point!). Thank you to [makeareadme.com](https://www.makeareadme.com/) for this template. - -## Suggestions for a good README -Every project is different, so consider which of these sections apply to yours. The sections used in the template are suggestions for most open source projects. Also keep in mind that while a README can be too long and detailed, too long is better than too short. If you think your README is too long, consider utilizing another form of documentation rather than cutting out information. - ## Name -Choose a self-explaining name for your project. +DAP is a Detector Analysis Pipeline which runs on detector stream data produced by sf-daq ## Description -Let people know what your project can do specifically. Provide context and add a link to any reference visitors might be unfamiliar with. A list of Features or a Background subsection can also be added here. If there are alternatives to your project, this is a good place to list differentiating factors. - -## Badges -On some READMEs, you may see small images that convey metadata, such as whether or not all the tests are passing for the project. You can use Shields to add some to your README. Many services also have instructions for adding a badge. - -## Visuals -Depending on what you are making, it can be a good idea to include screenshots or even a video (you'll frequently see GIFs rather than actual videos). Tools like ttygif can help, but check out Asciinema for a more sophisticated method. ## Installation -Within a particular ecosystem, there may be a common way of installing things, such as using Yarn, NuGet, or Homebrew. However, consider the possibility that whoever is reading your README is a novice and would like more guidance. Listing specific steps helps remove ambiguity and gets people to using your project as quickly as possible. If it only runs in a specific context like a particular programming language version or operating system or has dependencies that have to be installed manually, also add a Requirements subsection. +At PSI-Swissfel there is already pre-installed package/conda environment to use with: +``` +source /sf/jungfrau/applications/miniconda3/etc/profile.d/conda.sh +conda activate dap +``` + +### To install from source +Create conda environment +``` +conda create -n test-dap cython numpy pyzmq jungfrau_util +conda activate test-dap +``` +Clone code of the dap and install peakfinder8_extension in the conda environment +``` +git clone https://gitlab.psi.ch/sf-daq/dap.git +cd dap +make install +``` ## Usage Use examples liberally, and show the expected output if you can. It's helpful to have inline the smallest example of usage that you can demonstrate, while providing links to more sophisticated examples if they are too long to reasonably include in the README. -## Support -Tell people where they can go to for help. It can be any combination of an issue tracker, a chat room, an email address, etc. +## Acknowledgment -## Roadmap -If you have ideas for releases in the future, it is a good idea to list them in the README. - -## Contributing -State if you are open to contributions and what your requirements are for accepting them. - -For people who want to make changes to your project, it's helpful to have some documentation on how to get started. Perhaps there is a script that they should run or some environment variables that they need to set. Make these steps explicit. These instructions could also be useful to your future self. - -You can also document commands to lint the code or run tests. These steps help to ensure high code quality and reduce the likelihood that the changes inadvertently break something. Having instructions for running tests is especially helpful if it requires external setup, such as starting a Selenium server for testing in a browser. - -## Authors and acknowledgment -Show your appreciation to those who have contributed to the project. - -## License -For open source projects, say how it is licensed. - -## Project status -If you have run out of energy or time for your project, put a note at the top of the README saying that development has slowed down or stopped completely. Someone may choose to fork your project or volunteer to step in as a maintainer or owner, allowing your project to keep going. You can also make an explicit request for maintainers. + Initially DAP was made for SFX analysis, to run peak finder on a stream of data. peakfinder8 from cheetah was used for this purpose and big thanks to Valerio Mariani who provided cython implementation of peakfinder8 diff --git a/ap/accumulator.py b/ap/accumulator.py new file mode 100644 index 0000000..fdb7794 --- /dev/null +++ b/ap/accumulator.py @@ -0,0 +1,66 @@ +import zmq +import os +import argparse + +flags = 0 + +OUTPUT_DIR_NAME = "/gpfs/photonics/swissfel/buffer/dap/data" + +def main(): + + parser = argparse.ArgumentParser() + + parser.add_argument("--accumulator", default="localhost", help="name of host where accumulator works") + parser.add_argument("--accumulator_port", default=13002, type=int, help="accumulator port") + + args = parser.parse_args() + + FA_HOST_ACCUMULATE = args.accumulator + FA_PORT_ACCUMULATE = args.accumulator_port + + zmq_context = zmq.Context(io_threads=4) + poller = zmq.Poller() + +# Accumulator + if True: + accumulator_socket = zmq_context.socket(zmq.PULL) + accumulator_socket.bind('tcp://*:%s' % FA_PORT_ACCUMULATE ) + + poller.register(accumulator_socket, zmq.POLLIN) + + n_frames_received = 0 + + run_name_before = 'very_first_start' + fNameOutput = f'{OUTPUT_DIR_NAME}/{run_name_before}.dap' + if not os.path.isdir(os.path.dirname(fNameOutput)): + os.makedirs(os.path.dirname(fNameOutput)) + outputDap=open(fNameOutput, 'a') + + while True: + events = dict(poller.poll(10)) # in accumulator check for worker output every 0.01 seconds + + if accumulator_socket in events: + results = accumulator_socket.recv_json(flags) + n_frames_received += 1 + + pulse_id = results.get('pulse_id', 0) + run_name = str(pulse_id//10000*10000) + detector = results.get('detector_name', "") + + if run_name != run_name_before: + run_name_before = run_name + outputDap.close() + fNameOutput = f'{OUTPUT_DIR_NAME}/{detector}/{run_name_before}.dap' + if not os.path.isdir(os.path.dirname(fNameOutput)): + os.makedirs(os.path.dirname(fNameOutput)) + outputDap=open(fNameOutput, 'a') + + pr_rois = results.get("roi_intensities",[]) + print(pulse_id, results.get('is_good_frame', -1), results.get('is_hit_frame', False), results.get('number_of_spots', -1), results.get("laser_on", False), *pr_rois, file=outputDap) + + else: + outputDap.flush() # may be too intensive + +if __name__ == "__main__": + main() + diff --git a/ap/worker.py b/ap/worker.py new file mode 100644 index 0000000..80e2732 --- /dev/null +++ b/ap/worker.py @@ -0,0 +1,501 @@ +from copy import copy +import numpy as np +import zmq +from time import time, sleep +from datetime import datetime +import os +import argparse +from random import randint, gauss +from math import exp +import json + +from peakfinder8_extension import peakfinder_8 +import jungfrau_utils as ju + +flags = 0 + +MODULE_SIZE_X = 1024 +MODULE_SIZE_Y = 512 + +def radial_profile(data, r, nr, keep_pixels=None): + if keep_pixels is not None: + tbin = np.bincount(r, data[keep_pixels].ravel()) + else: + tbin = np.bincount(r, data.ravel()) + radialprofile = tbin / nr + return radialprofile + +def prepare_radial_profile(data, center, keep_pixels=None): + y, x = np.indices((data.shape)) + r = np.sqrt((x - center[0])**2 + (y - center[1])**2) + if keep_pixels is not None: + r = r[keep_pixels].astype(int).ravel() + else: + r = r.astype(np.int).ravel() + nr = np.bincount(r) + return r, nr + + +def mask_corner_pixels(pixel_mask, n_pixels): + + return pixel_mask + + +def main(): + + parser = argparse.ArgumentParser() + + parser.add_argument("--backend", default=None, help="backend address") + parser.add_argument("--accumulator", default="localhost", help="name of host where accumulator works") + parser.add_argument("--accumulator_port", default=13002, type=int, help="accumulator port") + parser.add_argument("--visualisation", default="localhost", help="name of host where visualisation works") + parser.add_argument("--visualisation_port", default=13002, type=int, help="visualisation port") + parser.add_argument("--peakfinder_parameters", default=None, help="json file with peakfinder parameters") + parser.add_argument("--skip_frames_rate", default=1, type=int, help="send to streamvis each of skip_frames_rate frames") + + args = parser.parse_args() + + if args.backend: + BACKEND_ADDRESS = args.backend + else: + print("no backend address defined, exit") + exit() + + FA_HOST_ACCUMULATE = args.accumulator + FA_PORT_ACCUMULATE = args.accumulator_port + FA_HOST_VISUALISATION = args.visualisation + FA_PORT_VISUALISATION = args.visualisation_port + + skip_frames_rate = args.skip_frames_rate + + peakfinder_parameters = {} + peakfinder_parameters_time = -1 + if args.peakfinder_parameters is not None and os.path.exists(args.peakfinder_parameters): + with open(args.peakfinder_parameters, "r") as read_file: + peakfinder_parameters = json.load(read_file) + peakfinder_parameters_time = os.path.getmtime(args.peakfinder_parameters) + + pulseid = 0 + + ju_stream_adapter = ju.StreamAdapter() + + zmq_context = zmq.Context(io_threads=4) + poller = zmq.Poller() + +# all the normal workers + if True: + worker = 1 + + # receive from backend: + backend_socket = zmq_context.socket(zmq.PULL) + backend_socket.connect(BACKEND_ADDRESS) + + poller.register(backend_socket, zmq.POLLIN) + + accumulator_socket = zmq_context.socket(zmq.PUSH) + accumulator_socket.connect('tcp://%s:%s' % (FA_HOST_ACCUMULATE, FA_PORT_ACCUMULATE) ) + + visualisation_socket = zmq_context.socket(zmq.PUB) + visualisation_socket.connect('tcp://%s:%s' % (FA_HOST_VISUALISATION, FA_PORT_VISUALISATION) ) + + # in case of problem with communication to visualisation, keep in 0mq buffer only few messages + visualisation_socket.set_hwm(10) + + keep_pixels = None + r_radial_integration = None + center_radial_integration = None + + results = {} + + pedestal_file_name_saved = None + + pixel_mask_corrected = None + pixel_mask_pf = None + + n_aggregated_images = 1 + data_summed = None + + while True: + +# check if peakfinder parameters changed and then re-read it + try: + if peakfinder_parameters_time > 0: + new_time = os.path.getmtime(args.peakfinder_parameters) + if ( new_time - peakfinder_parameters_time ) > 2.0: + old_peakfinder_parameters = peakfinder_parameters + sleep(0.5) + with open(args.peakfinder_parameters, "r") as read_file: + peakfinder_parameters = json.load(read_file) + peakfinder_parameters_time = new_time + center_radial_integration = None + if worker == 0: + print(f'({pulseid}) update peakfinder parameters {old_peakfinder_parameters}', flush=True) + print(f' --> {peakfinder_parameters}', flush=True) + print("",flush=True) + except: + print(f'({pulseid}) problem to read peakfinder parameters file, worker : {worker}', flush=True) + + events = dict(poller.poll(2000)) # check every 2 seconds in each worker + if backend_socket in events: + + metadata = backend_socket.recv_json(flags) + image = backend_socket.recv(flags, copy=False, track=False) + image = np.frombuffer(image, dtype=metadata['type']).reshape(metadata['shape']) + + results = copy(metadata) + if results['shape'][0] == 2 and results['shape'][1] == 2: + continue + + pulseid = results.get('pulse_id', 0) + results.update(peakfinder_parameters) + + detector = results.get('detector_name', "") + + results['laser_on'] = False + results['number_of_spots'] = 0 + results['is_hit_frame'] = False + + daq_rec = results.get("daq_rec",0) + event_laser = bool((daq_rec>>16)&1) + event_darkshot = bool((daq_rec>>17)&1) + event_fel = bool((daq_rec>>18)&1) + event_ppicker = bool((daq_rec>>19)&1) +# + if not event_darkshot: + results['laser_on'] = event_laser + +# special settings for Cristallina JF16 detector +# if detector in ["JF16T03V01"]: +# if not event_ppicker: +# continue + +# special settings for p20270, only few shots were opened by pulse-picker +# if detector in ["JF06T32V02"]: +# if event_ppicker: +# results['number_of_spots'] = 50 +# results['is_hit_frame'] = True + + double_pixels = results.get('double_pixels', "mask") + + pedestal_file_name = metadata.get("pedestal_name", None) + if pedestal_file_name is not None and pedestal_file_name != pedestal_file_name_saved: + n_corner_pixels_mask = results.get("n_corner_pixels_mask", 0) + pixel_mask_current = ju_stream_adapter.handler.pixel_mask + pixel_mask_current = mask_corner_pixels(pixel_mask_current, n_corner_pixels_mask) + ju_stream_adapter.handler.pixel_mask = pixel_mask_current + pedestal_file_name_saved = pedestal_file_name + + data = ju_stream_adapter.process(image, metadata, double_pixels=double_pixels) + + # 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 == "": + continue + + 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 + 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) + id_pixel_mask_2 = id(pixel_mask_corrected) + + if id_pixel_mask_1 != id_pixel_mask_2: + keep_pixels = None + r_radial_integration = None + if pixel_mask_corrected is not None: + #pixel_mask_corrected=np.ascontiguousarray(pixel_mask_corrected) + pixel_mask_pf = np.ascontiguousarray(pixel_mask_corrected).astype(np.int8, copy=False) + + else: + pixel_mask_pf = None +# + disabled_modules = results.get("disabled_modules", []) + +# add additional mask at the edge of modules for JF06T08 + apply_additional_mask = (results.get("apply_additional_mask", 0) == 1) + if detector == "JF06T08V04" and apply_additional_mask: + # edge pixels + pixel_mask_pf[67:1097,1063] = 0 + pixel_mask_pf[0:1030, 1100] = 0 + + pixel_mask_pf[1106:2136, 1131] = 0 + pixel_mask_pf[1039:2069, 1168] = 0 + + pixel_mask_pf[1039:2069, 1718] = 0 + pixel_mask_pf[1039:2069, 1681] = 0 + + pixel_mask_pf[1106:2136, 618] = 0 + + pixel_mask_pf[1106:2136, 581] = 0 + + pixel_mask_pf[67:1097,513] = 0 + + pixel_mask_pf[67:1097, 550] = 0 + + pixel_mask_pf[0:1030, 1650] = 0 + + pixel_mask_pf[0:1030, 1613] = 0 + + pixel_mask_pf[1106, 68:582] = 0 + + pixel_mask_pf[1096, 550:1064] = 0 + pixel_mask_pf[1106, 618:1132] = 0 + + pixel_mask_pf[1029, 1100:1614] = 0 + pixel_mask_pf[1039, 1168:1682] = 0 + + pixel_mask_pf[1039, 1718:2230] = 0 + + pixel_mask_pf[1096, 0:513] = 0 + + pixel_mask_pf[1029, 1650:2163] = 0 + + pixel_mask_pf[2068, 1168:2232] = 0 + + pixel_mask_pf[67,0:1063] = 0 + + #bad region in left bottom inner module + pixel_mask_pf[842:1097, 669:671] = 0 + + #second bad region in left bottom inner module + pixel_mask_pf[1094, 620:807] = 0 + + # vertical line in upper left bottom module + pixel_mask_pf[842:1072, 87:90] = 0 + + pixel_mask_pf[1794, 1503:1550] = 0 + + if detector == "JF17T16V01" and apply_additional_mask: + # mask module 11 + pixel_mask_pf[2619:3333,1577:2607] = 0 + + if pixel_mask_corrected is not None: + 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) + results['saturated_pixels'] = len(saturated_pixels_coordinates[0]) + results['saturated_pixels_x'] = saturated_pixels_coordinates[1].tolist() + results['saturated_pixels_y'] = saturated_pixels_coordinates[0].tolist() + +# pump probe analysis + do_radial_integration = results.get("do_radial_integration", 0) + + if (do_radial_integration != 0): + + data_copy_1 = np.copy(data) + + if keep_pixels is None and pixel_mask_pf is not None: + keep_pixels = pixel_mask_pf!=0 + if center_radial_integration is None: + center_radial_integration = [results['beam_center_x'], results['beam_center_y']] + r_radial_integration = 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_min_max = [int(np.min(r_radial_integration)), int(np.max(r_radial_integration))+1] + + + apply_threshold = results.get('apply_threshold', 0) + + if (apply_threshold != 0) and all ( k in results for k in ('threshold_min', 'threshold_max')): + threshold_min = float(results['threshold_min']) + threshold_max = float(results['threshold_max']) + data_copy_1[data_copy_1 < threshold_min] = np.nan + if threshold_max > threshold_min: + data_copy_1[data_copy_1 > threshold_max] = np.nan + + rp = radial_profile(data_copy_1, r_radial_integration, nr_radial_integration, keep_pixels) + + silent_region_min = results.get("radial_integration_silent_min", 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 + silent_region_max > silent_region_min and + 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]) + rp = rp/integral_silent_region + results['radint_normalised'] = [silent_region_min, silent_region_max] + + results['radint_I'] = list(rp[r_min_max[0]:]) + results['radint_q'] = r_min_max + + #copy image to work with peakfinder, just in case + d=np.copy(data) + + # make all masked pixels values nan's + if pixel_mask_pf is not None: + d[pixel_mask_pf != 1] = np.nan + + apply_threshold = results.get('apply_threshold', 0) + threshold_value_choice = results.get('threshold_value', "NaN") + threshold_value = 0 if threshold_value_choice == "0" else np.nan + if (apply_threshold != 0) and all ( k in results for k in ('threshold_min', 'threshold_max')): + threshold_min = float(results['threshold_min']) + threshold_max = float(results['threshold_max']) + d[d < threshold_min] = threshold_value + if threshold_max > threshold_min: + d[d > threshold_max] = threshold_value + + # if roi calculation request is present, make it + roi_x1 = results.get('roi_x1', []) + roi_x2 = results.get('roi_x2', []) + roi_y1 = results.get('roi_y1', []) + roi_y2 = results.get('roi_y2', []) + + if len(roi_x1) > 0 and len(roi_x1) == len(roi_x2) and len(roi_x1) == len(roi_y1) and len(roi_x1) == len(roi_y2): + roi_results = [0]*len(roi_x1) + roi_results_normalised = [0]*len(roi_x1) + + if pixel_mask_pf is not None: + + results['roi_intensities_x'] = [] + results['roi_intensities_proj_x'] = [] + + for iRoi in range(len(roi_x1)): + data_roi = np.copy(d[roi_y1[iRoi]:roi_y2[iRoi],roi_x1[iRoi]:roi_x2[iRoi]]) + + roi_results[iRoi] = np.nansum(data_roi) + if threshold_value_choice == "NaN": + roi_results_normalised[iRoi] = roi_results[iRoi]/((roi_y2[iRoi]-roi_y1[iRoi])*(roi_x2[iRoi]-roi_x1[iRoi])) + else: + roi_results_normalised[iRoi] = np.nanmean(data_roi) + + results['roi_intensities_x'].append([roi_x1[iRoi],roi_x2[iRoi]]) + results['roi_intensities_proj_x'].append(np.nansum(data_roi,axis=0).tolist()) + + results['roi_intensities'] = [ float(r) for r in roi_results] + results['roi_intensities_normalised'] = [ float(r) for r in roi_results_normalised ] + +# SPI analysis + do_spi_analysis = results.get("do_spi_analysis", 0) + + if (do_spi_analysis != 0) and 'roi_intensities_normalised' in results and len(results['roi_intensities_normalised']) >= 2: + + if 'spi_limit' in results and len(results['spi_limit']) == 2: + + number_of_spots = 0 + if results['roi_intensities_normalised'][0] >= results['spi_limit'][0]: + number_of_spots += 25 + if results['roi_intensities_normalised'][1] >= results['spi_limit'][1]: + number_of_spots += 50 + + results['number_of_spots'] = number_of_spots + if number_of_spots > 0: + results['is_hit_frame'] = True + +# in case all needed parameters are present, make peakfinding + do_peakfinder_analysis = results.get('do_peakfinder_analysis', 0) + if (do_peakfinder_analysis != 0) and pixel_mask_pf is not None and all ( k in results for k in ('beam_center_x', 'beam_center_y', 'hitfinder_min_snr', 'hitfinder_min_pix_count', 'hitfinder_adc_thresh')): + x_beam = results['beam_center_x'] - 0.5 # to coordinates where position of first pixel/point is 0.5, 0.5 + y_beam = results['beam_center_y'] - 0.5 # to coordinates where position of first pixel/point is 0.5, 0.5 + hitfinder_min_snr = results['hitfinder_min_snr'] + hitfinder_min_pix_count = int(results['hitfinder_min_pix_count']) + hitfinder_adc_thresh = results['hitfinder_adc_thresh'] + + asic_ny, asic_nx = d.shape + nasics_y, nasics_x = 1, 1 + hitfinder_max_pix_count = 100 + max_num_peaks = 10000 + + # usually don't need to change this value, rather robust + hitfinder_local_bg_radius= 20. + + #in case of further modification with the mask, make a new one, independent from real mask + maskPr = np.copy(pixel_mask_pf) + + y, x = np.indices(d.shape) + pix_r = np.sqrt((x-x_beam)**2 + (y-y_beam)**2) + + (peak_list_x, peak_list_y, peak_list_value) = peakfinder_8(max_num_peaks, + d.astype(np.float32), + maskPr.astype(np.int8), + pix_r.astype(np.float32), + asic_nx, asic_ny, + nasics_x, nasics_y, + hitfinder_adc_thresh, + hitfinder_min_snr, + hitfinder_min_pix_count, + hitfinder_max_pix_count, + hitfinder_local_bg_radius) + + + number_of_spots = len(peak_list_x) + results['number_of_spots'] = number_of_spots + if number_of_spots != 0: + results['spot_x'] = [-1.0] * number_of_spots + results['spot_y'] = [-1.0] * number_of_spots + results['spot_intensity'] = copy(peak_list_value) + for i in range(number_of_spots): + results['spot_x'][i] = peak_list_x[i]+0.5 + results['spot_y'][i] = peak_list_y[i]+0.5 + else: + results['spot_x'] = [] + results['spot_y'] = [] + results['spot_intensity'] = [] + + npeaks_threshold_hit = results.get('npeaks_threshold_hit', 15) + + if number_of_spots >= npeaks_threshold_hit: + results['is_hit_frame'] = True + + forceSendVisualisation = False + if data.dtype != np.uint16: + apply_threshold = results.get('apply_threshold', 0) + apply_aggregation = results.get('apply_aggregation', 0) + if apply_aggregation == 0: + data_summed = None + n_aggregated_images = 1 + if (apply_threshold != 0) or (apply_aggregation != 0 ): + if (apply_threshold != 0) and all ( k in results for k in ('threshold_min', 'threshold_max')): + threshold_min = float(results['threshold_min']) + threshold_max = float(results['threshold_max']) + data[data < threshold_min] = 0.0 + if threshold_max > threshold_min: + data[data > threshold_max] = 0.0 + if (apply_aggregation != 0 ) and 'aggregation_max' in results: + if data_summed is not None: + data += data_summed + n_aggregated_images += 1 + data_summed = data.copy() + data_summed[data == -np.nan] = -np.nan + results['aggregated_images'] = n_aggregated_images + results['worker'] = worker + if n_aggregated_images >= results['aggregation_max']: + forceSendVisualisation = True + data_summed = None + n_aggregated_images = 1 + data[pixel_mask_pf == 0] = np.NaN + + else: + data = image + + results['type'] = str(data.dtype) + results['shape'] = data.shape + + frame_number = metadata['frame'] +# + accumulator_socket.send_json(results, flags) + + if (apply_aggregation != 0 ) and 'aggregation_max' in results: + if forceSendVisualisation: + visualisation_socket.send_json(results, flags | zmq.SNDMORE) + visualisation_socket.send(data, flags, copy=True, track=True) + else: + data = np.empty((2,2), dtype=np.uint16) + results['type'] = str(data.dtype) + results['shape'] = data.shape + visualisation_socket.send_json(results, flags | zmq.SNDMORE) + visualisation_socket.send(data, flags, copy=True, track=True) + else: + if results['is_good_frame'] and (results['is_hit_frame'] or randint(1, skip_frames_rate) == 1): + visualisation_socket.send_json(results, flags | zmq.SNDMORE) + visualisation_socket.send(data, flags, copy=True, track=True) + else: + data = np.empty((2,2), dtype=np.uint16) + results['type'] = str(data.dtype) + results['shape'] = data.shape + visualisation_socket.send_json(results, flags | zmq.SNDMORE) + visualisation_socket.send(data, flags, copy=True, track=True) + + +if __name__ == "__main__": + main() + diff --git a/peakfinder8/.gitignore b/peakfinder8/.gitignore new file mode 100644 index 0000000..bf577bb --- /dev/null +++ b/peakfinder8/.gitignore @@ -0,0 +1,5 @@ +build/ +cython/peakfinder8_extension.cpp +peakfinder8_extension.cpython-*-x86_64-linux-gnu.so +libcheetah.so +peakfinders.o diff --git a/peakfinder8/cheetahmodules.h b/peakfinder8/cheetahmodules.h new file mode 100644 index 0000000..d2867aa --- /dev/null +++ b/peakfinder8/cheetahmodules.h @@ -0,0 +1,46 @@ +/* + * cheetahmodules.h + * cheetah + * + * Created by Anton Barty on 7/2/11. + * Copyright 2011 CFEL. All rights reserved. + * + */ + +int peakfinder3(tPeakList*, float*, char*, long, long, long, long, float, float, long, long, long); +int peakfinder6(tPeakList*, float*, char*, long, long, long, long, float, float, long, long, long, float); +int peakfinder8(tPeakList*, float*, char*, float*, long, long, long, long, float, float, long, long, long); +int peakfinder8old(tPeakList*, float*, char*, float*, long, long, long, long, float, float, long, long, long); +int killNearbyPeaks(tPeakList*, float ); + + +#include +// from detectorObject.h file +//-------------------------------------------------------------------------------------------------------------------- +/* + * Bits for pixel masks + * Oriented along conventions defined for CXI file format ( https://github.com/FilipeMaia/CXI/raw/master/cxi_file_format.pdf ) + * CONVENTIONS: + * - All options are dominantly inherited during assembly and pixel integration (see assembleImage.cpp) + * - The default value for all options is "false" + * */ +static const uint16_t PIXEL_IS_PERFECT = 0; // Remember to change this value if necessary after adding a new option +static const uint16_t PIXEL_IS_INVALID = 1; // bit 0 +static const uint16_t PIXEL_IS_SATURATED = 2; // bit 1 +static const uint16_t PIXEL_IS_HOT = 4; // bit 2 +static const uint16_t PIXEL_IS_DEAD = 8; // bit 3 +static const uint16_t PIXEL_IS_SHADOWED = 16; // bit 4 +static const uint16_t PIXEL_IS_IN_PEAKMASK = 32; // bit 5 +static const uint16_t PIXEL_IS_TO_BE_IGNORED = 64; // bit 6 +static const uint16_t PIXEL_IS_BAD = 128; // bit 7 +static const uint16_t PIXEL_IS_OUT_OF_RESOLUTION_LIMITS = 256; // bit 8 +static const uint16_t PIXEL_IS_MISSING = 512; // bit 9 +static const uint16_t PIXEL_IS_NOISY = 1024; // bit 10 +static const uint16_t PIXEL_IS_ARTIFACT_CORRECTED = 2048; // bit 11 +static const uint16_t PIXEL_FAILED_ARTIFACT_CORRECTION = 4096; // bit 12 +static const uint16_t PIXEL_IS_PEAK_FOR_HITFINDER = 8192; // bit 13 +static const uint16_t PIXEL_IS_PHOTON_BACKGROUND_CORRECTED = 16384; // bit 14 +static const uint16_t PIXEL_IS_ALL = PIXEL_IS_INVALID | PIXEL_IS_SATURATED | PIXEL_IS_HOT | PIXEL_IS_DEAD | PIXEL_IS_SHADOWED | PIXEL_IS_IN_PEAKMASK | PIXEL_IS_TO_BE_IGNORED | PIXEL_IS_BAD | PIXEL_IS_OUT_OF_RESOLUTION_LIMITS | PIXEL_IS_MISSING | PIXEL_IS_NOISY | PIXEL_IS_ARTIFACT_CORRECTED | PIXEL_FAILED_ARTIFACT_CORRECTION | PIXEL_IS_PEAK_FOR_HITFINDER | PIXEL_IS_PHOTON_BACKGROUND_CORRECTED; // all bits +// +inline bool isAnyOfBitOptionsSet(uint16_t value, uint16_t option) {return ((value & option)!=0);} +//-------------------------------------------------------------------------------------------------------------------- diff --git a/peakfinder8/cython/peakfinder8.hh b/peakfinder8/cython/peakfinder8.hh new file mode 100644 index 0000000..2d5a971 --- /dev/null +++ b/peakfinder8/cython/peakfinder8.hh @@ -0,0 +1 @@ +int peakfinder8(tPeakList*, float*, char*, float*, long, long, long, long, float, float, long, long, long); diff --git a/peakfinder8/cython/peakfinder8_extension.pyx b/peakfinder8/cython/peakfinder8_extension.pyx new file mode 100644 index 0000000..e9071d3 --- /dev/null +++ b/peakfinder8/cython/peakfinder8_extension.pyx @@ -0,0 +1,90 @@ +cimport numpy +from libcpp.vector cimport vector +from libc.stdlib cimport malloc, free + +# +# PEAKFINDER 8 +# + +cdef extern from "peakfinders.h": + + ctypedef struct tPeakList: + long nPeaks + long nHot + float peakResolution + float peakResolutionA + float peakDensity + float peakNpix + float peakTotal + int memoryAllocated + long nPeaks_max + + float *peak_maxintensity + float *peak_totalintensity + float *peak_sigma + float *peak_snr + float *peak_npix + float *peak_com_x + float *peak_com_y + long *peak_com_index + float *peak_com_x_assembled + float *peak_com_y_assembled + float *peak_com_r_assembled + float *peak_com_q + float *peak_com_res + + void allocatePeakList(tPeakList* peak_list, long max_num_peaks) + void freePeakList(tPeakList peak_list) + +cdef extern from "peakfinder8.hh": + + 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) + +def peakfinder_8(int max_num_peaks, + numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] data, + numpy.ndarray[numpy.int8_t, ndim=2, mode="c"] mask, + numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] pix_r, + long asic_nx, long asic_ny, + long nasics_x, long nasics_y, float adc_thresh, float hitfinder_min_snr, + long hitfinder_min_pix_count, long hitfinder_max_pix_count, + long hitfinder_local_bg_radius): + + cdef numpy.int8_t *mask_pointer = &mask[0,0] + cdef char *mask_char_pointer = mask_pointer + + cdef tPeakList peak_list + allocatePeakList(&peak_list, max_num_peaks) + + peakfinder8(&peak_list, &data[0,0], mask_char_pointer, &pix_r[0,0], + asic_nx, asic_ny, nasics_x, nasics_y, + adc_thresh, hitfinder_min_snr, hitfinder_min_pix_count, + hitfinder_max_pix_count, hitfinder_local_bg_radius) + + cdef int i + cdef float peak_x, peak_y, peak_value + cdef vector[double] peak_list_x + cdef vector[double] peak_list_y + cdef vector[double] peak_list_value + + num_peaks = peak_list.nPeaks + + if num_peaks > max_num_peaks: + num_peaks = max_num_peaks + + for i in range(0, num_peaks): + + peak_x = peak_list.peak_com_x[i] + peak_y = peak_list.peak_com_y[i] + peak_value = peak_list.peak_totalintensity[i] + + peak_list_x.push_back(peak_x) + peak_list_y.push_back(peak_y) + peak_list_value.push_back(peak_value) + + freePeakList(peak_list) + + return (peak_list_x, peak_list_y, peak_list_value) diff --git a/peakfinder8/peakfinders.cpp b/peakfinder8/peakfinders.cpp new file mode 100644 index 0000000..bddf238 --- /dev/null +++ b/peakfinder8/peakfinders.cpp @@ -0,0 +1,566 @@ +// +// peakfinders.cpp +// cheetah +// +// Created by Anton Barty on 23/3/13. +// +// + +#include +#include +#include +#include +#include + +#include "peakfinders.h" +#include "cheetahmodules.h" + + +/* + * Create arrays for remembering Bragg peak data + */ +void allocatePeakList(tPeakList *peak, long NpeaksMax) { + peak->nPeaks = 0; + peak->nPeaks_max = NpeaksMax; + peak->nHot = 0; + peak->peakResolution = 0; + peak->peakResolutionA = 0; + peak->peakDensity = 0; + peak->peakNpix = 0; + peak->peakTotal = 0; + + peak->peak_maxintensity = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_totalintensity = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_sigma = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_snr = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_npix = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_com_x = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_com_y = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_com_index = (long *) calloc(NpeaksMax, sizeof(long)); + peak->peak_com_x_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_q = (float *) calloc(NpeaksMax, sizeof(float)); + peak->peak_com_res = (float *) calloc(NpeaksMax, sizeof(float)); + peak->memoryAllocated = 1; +} + +/* + * Clean up Bragg peak arrays + */ +void freePeakList(tPeakList peak) { + free(peak.peak_maxintensity); + free(peak.peak_totalintensity); + free(peak.peak_sigma); + free(peak.peak_snr); + free(peak.peak_npix); + free(peak.peak_com_x); + free(peak.peak_com_y); + free(peak.peak_com_index); + free(peak.peak_com_x_assembled); + free(peak.peak_com_y_assembled); + free(peak.peak_com_r_assembled); + free(peak.peak_com_q); + free(peak.peak_com_res); + peak.memoryAllocated = 0; +} + + + + +/* + * Peakfinder 8 + * Version before modifications during Cherezov December 2014 LE80 + * Count peaks by searching for connected pixels above threshold + * 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) { + + // Derived values + long pix_nx = asic_nx*nasics_x; + long pix_ny = asic_ny*nasics_y; + long pix_nn = pix_nx*pix_ny; + //long asic_nn = asic_nx*asic_ny; + long hitfinderNpeaksMax = peaklist->nPeaks_max; + + + peaklist->nPeaks = 0; + peaklist->peakNpix = 0; + peaklist->peakTotal = 0; + + + // Variables for this hitfinder + long nat = 0; + long lastnat = 0; + //long counter=0; + float total; + int search_x[] = {0,-1,0,1,-1,1,-1,0,1}; + int search_y[] = {0,-1,-1,-1,0,0,1,1,1}; + int search_n = 9; + long e; + long *inx = (long *) calloc(pix_nn, sizeof(long)); + long *iny = (long *) calloc(pix_nn, sizeof(long)); + float thisI, thisIraw; + float totI,totIraw; + float maxI, maxIraw; + float snr; + float peak_com_x; + float peak_com_y; + long thisx; + long thisy; + long fs, ss; + float com_x, com_y, com_e; + float thisADCthresh; + + + nat = 0; + //counter = 0; + total = 0.0; + snr=0; + maxI = 0; + + /* + * Create a buffer for image data so we don't nuke the main image by mistake + */ + float *temp = (float*) malloc(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) + */ + for(long i=0;i fmaxr) + fmaxr = pix_r[i]; + if (pix_r[i] < fminr) + fminr = pix_r[i]; + } + lmaxr = (int)ceil(fmaxr)+1; + lminr = 0; + + // Allocate and zero arrays + float *rsigma = (float*) malloc(lmaxr*sizeof(float)); + float *roffset = (float*) malloc(lmaxr*sizeof(float)); + long *rcount = (long*) malloc(lmaxr*sizeof(long)); + float *rthreshold = (float*) malloc(lmaxr*sizeof(float)); + + long *peakpixels = (long *) calloc(hitfinderMaxPixCount, sizeof(long)); + char *peakpixel = (char *) calloc(pix_nn, sizeof(char)); + + char *rthreshold_changed = (char *) malloc(lmaxr*sizeof(char)); + + int *pix_rint = (int *) malloc(pix_nn*sizeof(int)); + long *pixels_check = (long *) malloc(pix_nn*sizeof(long)); + + long peakCounter = 0; + + for(long i=0; i 0.1*rsigma[i] ) { + rthreshold_changed[i] = 1; + rthreshold_converged = false; + } + rthreshold[i] = thisthreshold; + } + } + } + + com_x=0; + com_y=0; + + //goto END; + + for(long mj=0; mj pix_nn) { + printf("Array bounds error: e=%li\n",e); + exit(1); + } + + thisr = pix_rint[e]; + thisADCthresh = rthreshold[thisr]; + + if(temp[e] > thisADCthresh && peakpixel[e] == 0){ + // This might be the start of a new peak - start searching + inx[0] = i; + iny[0] = j; + peakpixels[0] = e; + nat = 1; + totI = 0; + totIraw = 0; + maxI = 0; + maxIraw = 0; + peak_com_x = 0; + peak_com_y = 0; + + // Keep looping until the pixel count within this peak does not change + do { + + lastnat = nat; + // Loop through points known to be within this peak + for(long p=0; p= asic_nx) + continue; + if((iny[p]+search_y[k]) < 0) + continue; + if((iny[p]+search_y[k]) >= asic_ny) + continue; + + // Neighbour point in big array + thisx = inx[p]+search_x[k]+mi*asic_nx; + thisy = iny[p]+search_y[k]+mj*asic_ny; + e = thisx + thisy*pix_nx; + + //if(e < 0 || e >= pix_nn){ + // printf("Array bounds error: e=%i\n",e); + // continue; + //} + + thisr = pix_rint[e]; + thisADCthresh = rthreshold[thisr]; + + // Above threshold? + if(temp[e] > thisADCthresh && peakpixel[e] == 0 && mask[e] != 0){ + //if(nat < 0 || nat >= global->pix_nn) { + // printf("Array bounds error: nat=%i\n",nat); + // break + //} + thisI = temp[e] - roffset[thisr]; + totI += thisI; // add to integrated intensity + totIraw += temp[e]; + peak_com_x += thisI*( (float) thisx ); // for center of mass x + peak_com_y += thisI*( (float) thisy ); // for center of mass y + //temp[e] = 0; // zero out this intensity so that we don't count it again + inx[nat] = inx[p]+search_x[k]; + iny[nat] = iny[p]+search_y[k]; + peakpixel[e] = 1; + if(nat < hitfinderMaxPixCount) + peakpixels[nat] = e; + if (thisI > maxI) + maxI = thisI; + if (thisI > maxIraw) + maxIraw = temp[e]; + + nat++; + } + } + } + } while(lastnat != nat); + + + // Too many or too few pixels means ignore this 'peak'; move on now + if(nathitfinderMaxPixCount) { + continue; + } + + + /* + * Calculate center of mass for this peak from initial peak search + */ + com_x = peak_com_x/fabs(totI); + com_y = peak_com_y/fabs(totI); + com_e = lrint(com_x) + lrint(com_y)*pix_nx; + + long com_xi = lrint(com_x) - mi*asic_nx; + 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 + * (excluding pixels which look like they might be part of another peak) + */ + float localSigma=0; + float localOffset=0; + long ringWidth = 2*hitfinderLocalBGRadius; + + float sumI = 0; + float sumIsquared = 0; + long np_sigma = 0; + long np_counted = 0; + float fbgr; + float backgroundMaxI=0; + float fBackgroundThresh=0; + + for(long bj=-ringWidth; bj= asic_nx) + continue; + if((com_yi+bj) < 0) + continue; + if((com_yi+bj) >= asic_ny) + continue; + + // Within outer ring check + fbgr = sqrt( bi*bi + bj*bj ); + if( fbgr > ringWidth )// || fbgr <= hitfinderLocalBGRadius ) // || fbgr > hitfinderLocalBGRadius) + continue; + + // Position of this point in data stream + thisx = com_xi + bi + mi*asic_nx; + thisy = com_yi + bj + mj*asic_ny; + e = thisx + thisy*pix_nx; + + thisr = pix_rint[e]; + thisADCthresh = rthreshold[thisr]; + + // Intensity above background + thisI = temp[e]; + + + // If above ADC threshold, this could be part of another peak + //if (temp[e] > thisADCthresh) + // continue; + + // Keep track of value and value-squared for offset and sigma calculation + // if(peakpixel[e] == 0 && mask[e] != 0) { + if(temp[e] < thisADCthresh && peakpixel[e] == 0 && mask[e] != 0) { + np_sigma++; + sumI += thisI; + sumIsquared += (thisI*thisI); + if(thisI > backgroundMaxI) { + backgroundMaxI = thisI; + } + } + np_counted += 1; + } + } + + // Calculate local background and standard deviation + if (np_sigma != 0) { + localOffset = sumI/np_sigma; + localSigma = sqrt(sumIsquared/np_sigma - ((sumI/np_sigma)*(sumI/np_sigma))); + } + else { + localOffset = roffset[pix_rint[lrint(com_e)]]; + localSigma = 0.01; + } + + + /* + * Re-integrate (and re-centroid) peak using local background estimates + */ + totI = 0; + totIraw = 0; + maxI = 0; + maxIraw = 0; + peak_com_x = 0; + peak_com_y = 0; + for(long counter=1; counter maxIraw) + maxIraw = thisIraw; + if (thisI > maxI) + maxI = thisI; + } + com_x = peak_com_x/fabs(totI); + com_y = peak_com_y/fabs(totI); + com_e = lrint(com_x) + lrint(com_y)*pix_nx; + + + + /* + * Calculate signal-to-noise and apply SNR criteria + */ + snr = (float) (totI)/localSigma; + //snr = (float) (maxI)/localSigma; + //snr = (float) (totIraw-nat*localOffset)/localSigma; + //snr = (float) (maxIraw-localOffset)/localSigma; + + // The more pixels there are in the peak, the more relaxed we are about this criterion + if( snr < hitfinderMinSNR ) // - nat +hitfinderMinPixCount + continue; + + // 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 + //fBackgroundThresh = hitfinderMinSNR - nat; + //if(fBackgroundThresh > 4) fBackgroundThresh = 4; + fBackgroundThresh = 1; + fBackgroundThresh *= (backgroundMaxI-localOffset); + if( maxI < fBackgroundThresh) + continue; + + + // This is a peak? If so, add info to peak list + if(nat>=hitfinderMinPixCount && nat<=hitfinderMaxPixCount ) { + + // This CAN happen! + if(totI == 0) + continue; + + //com_x = peak_com_x/fabs(totI); + //com_y = peak_com_y/fabs(totI); + + e = lrint(com_x) + lrint(com_y)*pix_nx; + if(e < 0 || e >= pix_nn){ + printf("Array bounds error: e=%ld\n",e); + continue; + } + + // Remember peak information + if (peakCounter < hitfinderNpeaksMax) { + peaklist->peakNpix += nat; + peaklist->peakTotal += totI; + peaklist->peak_com_index[peakCounter] = e; + peaklist->peak_npix[peakCounter] = nat; + peaklist->peak_com_x[peakCounter] = com_x; + peaklist->peak_com_y[peakCounter] = com_y; + peaklist->peak_totalintensity[peakCounter] = totI; + peaklist->peak_maxintensity[peakCounter] = maxI; + peaklist->peak_sigma[peakCounter] = localSigma; + peaklist->peak_snr[peakCounter] = snr; + peakCounter++; + peaklist->nPeaks = peakCounter; + } + else { + peakCounter++; + } + } + } + } + } + } + } + + //END: ; + + free(temp); + free(inx); + free(iny); + free(peakpixel); + free(peakpixels); + + + free(roffset); + free(rsigma); + free(rcount); + free(rthreshold); + + free(pix_rint); + free(pixels_check); + free(rthreshold_changed); + + peaklist->nPeaks = peakCounter; + return(peaklist->nPeaks); + /*************************************************/ + + +} + + diff --git a/peakfinder8/peakfinders.h b/peakfinder8/peakfinders.h new file mode 100644 index 0000000..eac0bf8 --- /dev/null +++ b/peakfinder8/peakfinders.h @@ -0,0 +1,54 @@ +// +// peakfinders.h +// cheetah +// +// Created by Anton Barty on 23/3/13. +// +// + +#ifndef cheetah_peakfinders_h +#define cheetah_peakfinders_h + + + +typedef struct { +public: + long nPeaks; + long nHot; + float peakResolution; // Radius of 80% of peaks + float peakResolutionA; // Radius of 80% of peaks + float peakDensity; // Density of peaks within this 80% figure + float peakNpix; // Number of pixels in peaks + float peakTotal; // Total integrated intensity in peaks + int memoryAllocated; + long nPeaks_max; + + float *peak_maxintensity; // Maximum intensity in peak + float *peak_totalintensity; // Integrated intensity in peak + float *peak_sigma; // Signal-to-noise ratio of peak + float *peak_snr; // Signal-to-noise ratio of peak + float *peak_npix; // Number of pixels in peak + float *peak_com_x; // peak center of mass x (in raw layout) + float *peak_com_y; // peak center of mass y (in raw layout) + long *peak_com_index; // closest pixel corresponding to peak + float *peak_com_x_assembled; // peak center of mass x (in assembled layout) + float *peak_com_y_assembled; // peak center of mass y (in assembled layout) + float *peak_com_r_assembled; // peak center of mass r (in assembled layout) + float *peak_com_q; // Scattering vector of this peak + float *peak_com_res; // REsolution of this peak +} tPeakList; + + + +//void cPeakList::cPeakList(void) { +//} + +//void cPeakList::~cPeakList(void) { +// free(); +//} + +void allocatePeakList(tPeakList*, long); +void freePeakList(tPeakList); + + +#endif diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..90039b3 --- /dev/null +++ b/setup.py @@ -0,0 +1,15 @@ +peakfinder8_include_dir = 'peakfinder8/' +peakfinder8_library_dir = 'peakfinder8/' + +from distutils.core import setup , Extension +from Cython.Build import cythonize +import numpy + +peakfinder8_ext = Extension ( "peakfinder8_extension" , sources = [ "peakfinder8/cython/peakfinder8_extension.pyx"] , + include_dirs = [ peakfinder8_include_dir, numpy.get_include() ], + library_dirs = [ peakfinder8_library_dir ], + libraries=["peakfinder8"], + language = "c++" ) + + +setup ( name = "peakfinder8_extension" , ext_modules = cythonize ( peakfinder8_ext ))