diff --git a/Makefile b/Makefile index e1e9aa2..e608028 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ clean: 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 +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 diff --git a/README.md b/README.md index e76770d..794bc42 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ Clone and install dap git clone https://gitlab.psi.ch/sf-daq/dap.git cd dap make install -``` +``` # Architecture @@ -37,7 +37,7 @@ The dap architecture is designed for horizontal scalability, processing various - Metadata-only results to the accumulator for storage. ## 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: ```bash python dap/worker.py --help @@ -65,7 +65,7 @@ options: 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. ## Accumulator @@ -87,9 +87,9 @@ options: ``` # 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. 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. * **Radial Profile Integration** - + 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. * `'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. * `'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_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_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_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** @@ -156,10 +156,10 @@ options: * 0: if intensity in both ROIs falls below the respective thresholds * 25: ROI1 has high energy but not ROI2 * 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** - + 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: @@ -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. Algorithm Output: - * `'aggregated_images': int` - Indicates the count of aggregated images. - + * `'aggregated_images': int` - Indicates the count of aggregated images. + * **Frame Tagging** When event propagation is integrated into the detector, dap allows frames to be tagged accordingly, facilitating their categorization during visualization. - + 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". - + * **Detection of Saturated Pixels** For every frame received by dap, an analysis is performed to ascertain the quantity and positions of saturated pixels. - + Algorithm Output: * `'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** @@ -194,12 +194,12 @@ options: * **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. - + Use the `'apply_additional_mask': 0/1` - Input flag to enable this functionality. * **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. ## Input parameters (File) diff --git a/dap/accumulator.py b/dap/accumulator.py index fdb7794..ed65115 100644 --- a/dap/accumulator.py +++ b/dap/accumulator.py @@ -14,14 +14,14 @@ def main(): 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 +# Accumulator if True: accumulator_socket = zmq_context.socket(zmq.PULL) accumulator_socket.bind('tcp://*:%s' % FA_PORT_ACCUMULATE ) @@ -42,7 +42,7 @@ def main(): 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', "") diff --git a/dap/worker.py b/dap/worker.py index 6c1ad34..d6ba37e 100644 --- a/dap/worker.py +++ b/dap/worker.py @@ -23,7 +23,7 @@ def radial_profile(data, r, nr, keep_pixels=None): else: tbin = np.bincount(r, data.ravel()) radialprofile = tbin / nr - return radialprofile + return radialprofile def prepare_radial_profile(data, center, keep_pixels=None): 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") args = parser.parse_args() - + if args.backend: BACKEND_ADDRESS = args.backend else: @@ -87,7 +87,7 @@ def main(): worker = 1 # receive from backend: - backend_socket = zmq_context.socket(zmq.PULL) + backend_socket = zmq_context.socket(zmq.PULL) backend_socket.connect(BACKEND_ADDRESS) poller.register(backend_socket, zmq.POLLIN) @@ -140,7 +140,7 @@ def main(): 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']) + image = np.frombuffer(image, dtype=metadata['type']).reshape(metadata['shape']) results = copy(metadata) if results['shape'][0] == 2 and results['shape'][1] == 2: @@ -160,7 +160,7 @@ def main(): 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 @@ -174,7 +174,7 @@ def main(): # 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) @@ -190,12 +190,12 @@ def main(): # 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 + # 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) + 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: @@ -206,14 +206,14 @@ def main(): pixel_mask_pf = np.ascontiguousarray(pixel_mask_corrected).astype(np.int8, copy=False) else: - pixel_mask_pf = None + 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 + # edge pixels pixel_mask_pf[67:1097,1063] = 0 pixel_mask_pf[0:1030, 1100] = 0 @@ -232,7 +232,7 @@ def main(): 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 @@ -267,7 +267,7 @@ def main(): 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) @@ -277,7 +277,7 @@ def main(): # pump probe analysis do_radial_integration = results.get("do_radial_integration", 0) - + if (do_radial_integration != 0): data_copy_1 = np.copy(data) @@ -290,9 +290,9 @@ def main(): 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) + + 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']) @@ -305,11 +305,11 @@ def main(): 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 + + 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] @@ -352,7 +352,7 @@ def main(): 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) + 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: @@ -377,7 +377,7 @@ def main(): if results['roi_intensities_normalised'][1] >= results['spi_limit'][1]: number_of_spots += 50 - results['number_of_spots'] = number_of_spots + results['number_of_spots'] = number_of_spots if number_of_spots > 0: results['is_hit_frame'] = True diff --git a/peakfinder8/cheetahmodules.h b/peakfinder8/cheetahmodules.h index d2867aa..ef28cda 100644 --- a/peakfinder8/cheetahmodules.h +++ b/peakfinder8/cheetahmodules.h @@ -16,7 +16,7 @@ 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 ) diff --git a/peakfinder8/peakfinders.cpp b/peakfinder8/peakfinders.cpp index bddf238..d6b64ec 100644 --- a/peakfinder8/peakfinders.cpp +++ b/peakfinder8/peakfinders.cpp @@ -21,14 +21,14 @@ */ void allocatePeakList(tPeakList *peak, long NpeaksMax) { peak->nPeaks = 0; - peak->nPeaks_max = NpeaksMax; + 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)); @@ -41,7 +41,7 @@ void allocatePeakList(tPeakList *peak, long NpeaksMax) { 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->peak_com_res = (float *) calloc(NpeaksMax, sizeof(float)); peak->memoryAllocated = 1; } @@ -75,20 +75,20 @@ void freePeakList(tPeakList peak) { * 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; @@ -111,28 +111,28 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long 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) @@ -150,16 +150,16 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long } 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)); @@ -177,7 +177,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long pixels_check[i] = i; } long n_pixels_check = pix_nn; - + // Compute sigma and average of data values 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) @@ -215,7 +215,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long rcount[thisr] += 1; } 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_y=0; @@ -257,24 +257,24 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long 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; @@ -287,10 +287,10 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long 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_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) { @@ -340,20 +340,20 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long 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 */ @@ -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_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) @@ -372,7 +372,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long float localSigma=0; float localOffset=0; long ringWidth = 2*hitfinderLocalBGRadius; - + float sumI = 0; float sumIsquared = 0; long np_sigma = 0; @@ -380,10 +380,10 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long float fbgr; float backgroundMaxI=0; float fBackgroundThresh=0; - + for(long bj=-ringWidth; 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) { @@ -428,7 +428,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long np_counted += 1; } } - + // Calculate local background and standard deviation if (np_sigma != 0) { 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)]]; localSigma = 0.01; } - - + + /* * 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]; thisIraw = temp[e]; thisI = thisIraw - localOffset; - + totI += thisI; 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_e = lrint(com_x) + lrint(com_y)*pix_nx; - + /* * 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) (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; @@ -495,24 +495,24 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long 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; @@ -539,28 +539,28 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long } //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 index eac0bf8..7d3a371 100644 --- a/peakfinder8/peakfinders.h +++ b/peakfinder8/peakfinders.h @@ -22,7 +22,7 @@ public: 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