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

more consistent whitespace/aligment

Browse Source
This commit is contained in:
augustin_s
2024-03-22 17:29:35 +01:00
parent d2ce6a8e16
commit 80620d1f14
4 changed files with 198 additions and 173 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
peakfinder8/cheetahmodules.h
View File

@ -24,23 +24,23 @@ int killNearbyPeaks(tPeakList*, float );
* - All options are dominantly inherited during assembly and pixel integration (see assembleImage.cpp) * - All options are dominantly inherited during assembly and pixel integration (see assembleImage.cpp)
* - The default value for all options is "false" * - 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_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_INVALID = 1; // bit 0
static const uint16_t PIXEL_IS_SATURATED = 2; // bit 1 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_HOT = 4; // bit 2
static const uint16_t PIXEL_IS_DEAD = 8; // bit 3 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_SHADOWED = 16; // bit 4
static const uint16_t PIXEL_IS_IN_PEAKMASK = 32; // bit 5 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_TO_BE_IGNORED = 64; // bit 6
static const uint16_t PIXEL_IS_BAD = 128; // bit 7 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_OUT_OF_RESOLUTION_LIMITS = 256; // bit 8
static const uint16_t PIXEL_IS_MISSING = 512; // bit 9 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_NOISY = 1024; // bit 10
static const uint16_t PIXEL_IS_ARTIFACT_CORRECTED = 2048; // bit 11 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_FAILED_ARTIFACT_CORRECTION = 4096; // bit 12
static const uint16_t PIXEL_IS_PEAK_FOR_HITFINDER = 8192; // bit 13 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_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 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);} inline bool isAnyOfBitOptionsSet(uint16_t value, uint16_t option) {return ((value & option)!=0);}
//-------------------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------------------

106
peakfinder8/cython/peakfinder8_extension.pyx
View File

@ -9,49 +9,66 @@ from libc.stdlib cimport malloc, free
cdef extern from "peakfinders.h": cdef extern from "peakfinders.h":
ctypedef struct tPeakList: ctypedef struct tPeakList:
long nPeaks long nPeaks
long nHot long nHot
float peakResolution float peakResolution
float peakResolutionA float peakResolutionA
float peakDensity float peakDensity
float peakNpix float peakNpix
float peakTotal float peakTotal
int memoryAllocated int memoryAllocated
long nPeaks_max long nPeaks_max
float *peak_maxintensity float *peak_maxintensity
float *peak_totalintensity float *peak_totalintensity
float *peak_sigma float *peak_sigma
float *peak_snr float *peak_snr
float *peak_npix float *peak_npix
float *peak_com_x float *peak_com_x
float *peak_com_y float *peak_com_y
long *peak_com_index long *peak_com_index
float *peak_com_x_assembled float *peak_com_x_assembled
float *peak_com_y_assembled float *peak_com_y_assembled
float *peak_com_r_assembled float *peak_com_r_assembled
float *peak_com_q float *peak_com_q
float *peak_com_res float *peak_com_res
void allocatePeakList(tPeakList* peak_list, long max_num_peaks) void allocatePeakList(tPeakList* peak_list, long max_num_peaks)
void freePeakList(tPeakList peak_list) void freePeakList(tPeakList peak_list)
cdef extern from "peakfinder8.hh": cdef extern from "peakfinder8.hh":
int peakfinder8(tPeakList *peaklist, int peakfinder8(
float *data, char *mask, float *pix_r, long asic_nx, long asic_ny, tPeakList *peaklist,
long nasics_x, long nasics_y, float ADCthresh, float hitfinderMinSNR, float *data,
long hitfinderMinPixCount, long hitfinderMaxPixCount, char *mask,
long hitfinderLocalBGRadius) 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, def peakfinder_8(
numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] data, int max_num_peaks,
numpy.ndarray[numpy.int8_t, ndim=2, mode="c"] mask, numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] data,
numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] pix_r, numpy.ndarray[numpy.int8_t, ndim=2, mode="c"] mask,
long asic_nx, long asic_ny, numpy.ndarray[numpy.float32_t, ndim=2, mode="c"] pix_r,
long nasics_x, long nasics_y, float adc_thresh, float hitfinder_min_snr, long asic_nx,
long hitfinder_min_pix_count, long hitfinder_max_pix_count, long asic_ny,
long hitfinder_local_bg_radius): 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 numpy.int8_t *mask_pointer = &mask[0,0]
cdef char *mask_char_pointer = <char*> mask_pointer cdef char *mask_char_pointer = <char*> mask_pointer
@ -59,10 +76,21 @@ def peakfinder_8(int max_num_peaks,
cdef tPeakList peak_list cdef tPeakList peak_list
allocatePeakList(&peak_list, max_num_peaks) allocatePeakList(&peak_list, max_num_peaks)
peakfinder8(&peak_list, &data[0,0], mask_char_pointer, &pix_r[0,0], peakfinder8(
asic_nx, asic_ny, nasics_x, nasics_y, &peak_list,
adc_thresh, hitfinder_min_snr, hitfinder_min_pix_count, &data[0,0],
hitfinder_max_pix_count, hitfinder_local_bg_radius) 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 int i
cdef float peak_x, peak_y, peak_value cdef float peak_x, peak_y, peak_value

223
peakfinder8/peakfinders.cpp
View File

@ -20,54 +20,53 @@
* Create arrays for remembering Bragg peak data * Create arrays for remembering Bragg peak data
*/ */
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));
peak->peak_snr = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_snr = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_npix = (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_x = (float *) calloc(NpeaksMax, sizeof(float));
peak->peak_com_y = (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_index = (long *) calloc(NpeaksMax, sizeof(long));
peak->peak_com_x_assembled = (float *) calloc(NpeaksMax, sizeof(float)); peak->peak_com_x_assembled = (float *) calloc(NpeaksMax, sizeof(float));
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;
} }
/* /*
* Clean up Bragg peak arrays * Clean up Bragg peak arrays
*/ */
void freePeakList(tPeakList peak) { void freePeakList(tPeakList peak) {
free(peak.peak_maxintensity); free(peak.peak_maxintensity);
free(peak.peak_totalintensity); free(peak.peak_totalintensity);
free(peak.peak_sigma); free(peak.peak_sigma);
free(peak.peak_snr); free(peak.peak_snr);
free(peak.peak_npix); free(peak.peak_npix);
free(peak.peak_com_x); free(peak.peak_com_x);
free(peak.peak_com_y); free(peak.peak_com_y);
free(peak.peak_com_index); free(peak.peak_com_index);
free(peak.peak_com_x_assembled); free(peak.peak_com_x_assembled);
free(peak.peak_com_y_assembled); free(peak.peak_com_y_assembled);
free(peak.peak_com_r_assembled); free(peak.peak_com_r_assembled);
free(peak.peak_com_q); free(peak.peak_com_q);
free(peak.peak_com_res); free(peak.peak_com_res);
peak.memoryAllocated = 0; peak.memoryAllocated = 0;
} }
/* /*
* Peakfinder 8 * Peakfinder 8
* Version before modifications during Cherezov December 2014 LE80 * Version before modifications during Cherezov December 2014 LE80
@ -100,7 +99,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
long e; long e;
long *inx = (long *) calloc(pix_nn, sizeof(long)); long *inx = (long *) calloc(pix_nn, sizeof(long));
long *iny = (long *) calloc(pix_nn, sizeof(long)); long *iny = (long *) calloc(pix_nn, sizeof(long));
float thisI, thisIraw; float thisI, thisIraw;
float totI,totIraw; float totI,totIraw;
float maxI, maxIraw; float maxI, maxIraw;
float snr; float snr;
@ -160,100 +159,100 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
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));
long *pixels_check = (long *) malloc(pix_nn*sizeof(long)); long *pixels_check = (long *) malloc(pix_nn*sizeof(long));
long peakCounter = 0; long peakCounter = 0;
for(long i=0; i<lmaxr; i++) { for(long i=0; i<lmaxr; i++) {
rthreshold[i] = 1e9; rthreshold[i] = 1e9;
rthreshold_changed[i] = 1; rthreshold_changed[i] = 1;
} }
for(long i=0;i<pix_nn;i++){ for(long i=0;i<pix_nn;i++){
pix_rint[i] = lrint(pix_r[i]); pix_rint[i] = lrint(pix_r[i]);
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)
long thisr; long thisr;
float thisoffset, thissigma; float thisoffset, thissigma;
float thisthreshold; float thisthreshold;
int counter = 0; int counter = 0;
bool rthreshold_converged = false; bool rthreshold_converged = false;
//goto END; //goto END;
// for(long counter=0; counter<5; counter++) { // for(long counter=0; counter<5; counter++) {
while ( !rthreshold_converged & counter < 10 ) { while ( !rthreshold_converged & counter < 10 ) {
//printf("counter %i %i\n", counter, n_pixels_check); //printf("counter %i %i\n", counter, n_pixels_check);
counter++; counter++;
//for(long i=0; i<lmaxr; i++) { //for(long i=0; i<lmaxr; i++) {
// roffset[i] = 0; // roffset[i] = 0;
// rsigma[i] = 0; // rsigma[i] = 0;
// rcount[i] = 0; // rcount[i] = 0;
//} //}
memset(roffset,0,lmaxr*sizeof(float)); memset(roffset,0,lmaxr*sizeof(float));
memset(rsigma, 0,lmaxr*sizeof(float)); memset(rsigma, 0,lmaxr*sizeof(float));
memset(rcount, 0,lmaxr*sizeof(long)); memset(rcount, 0,lmaxr*sizeof(long));
long new_pixels_check=0; long new_pixels_check=0;
//for(long i=0;i<pix_nn;i++){ //for(long i=0;i<pix_nn;i++){
for(long i_pixel=0; i_pixel<n_pixels_check; i_pixel++) { for(long i_pixel=0; i_pixel<n_pixels_check; i_pixel++) {
long i = pixels_check[i_pixel]; long i = pixels_check[i_pixel];
thisr = pix_rint[i]; thisr = pix_rint[i];
if ( rthreshold_changed[thisr] == 1 ) { if ( rthreshold_changed[thisr] == 1 ) {
if(mask[i] != 0) { if(mask[i] != 0) {
if(temp[i] < rthreshold[thisr]) { if(temp[i] < rthreshold[thisr]) {
roffset[thisr] += temp[i]; roffset[thisr] += temp[i];
rsigma[thisr] += (temp[i]*temp[i]); rsigma[thisr] += (temp[i]*temp[i]);
rcount[thisr] += 1; rcount[thisr] += 1;
} }
pixels_check[new_pixels_check] = i; pixels_check[new_pixels_check] = i;
new_pixels_check++; new_pixels_check++;
} }
} }
} }
n_pixels_check = new_pixels_check; n_pixels_check = new_pixels_check;
rthreshold_converged = true; rthreshold_converged = true;
for(long i=0; i<lmaxr; i++) { for(long i=0; i<lmaxr; i++) {
if ( rthreshold_changed[i] == 1 ) { if ( rthreshold_changed[i] == 1 ) {
if(rcount[i] == 0) { if(rcount[i] == 0) {
roffset[i] = 0; roffset[i] = 0;
rsigma[i] = 0; rsigma[i] = 0;
thisthreshold = 1e9; thisthreshold = 1e9;
//rthreshold[i] = ADCthresh; // For testing //rthreshold[i] = ADCthresh; // For testing
} }
else { else {
thisoffset = roffset[i]/rcount[i]; thisoffset = roffset[i]/rcount[i];
thissigma = sqrt(rsigma[i]/rcount[i] - (thisoffset)*(thisoffset)); thissigma = sqrt(rsigma[i]/rcount[i] - (thisoffset)*(thisoffset));
roffset[i] = thisoffset; roffset[i] = thisoffset;
rsigma[i] = thissigma; rsigma[i] = thissigma;
thisthreshold = roffset[i] + hitfinderMinSNR*rsigma[i]; thisthreshold = roffset[i] + hitfinderMinSNR*rsigma[i];
if(thisthreshold < ADCthresh) if(thisthreshold < ADCthresh)
thisthreshold = ADCthresh; thisthreshold = ADCthresh;
//rthreshold[i] = ADCthresh; // For testing //rthreshold[i] = ADCthresh; // For testing
} }
rthreshold_changed[i] = 0; rthreshold_changed[i] = 0;
if ( fabs(thisthreshold-rthreshold[i]) > 0.1*rsigma[i] ) { if ( fabs(thisthreshold-rthreshold[i]) > 0.1*rsigma[i] ) {
rthreshold_changed[i] = 1; rthreshold_changed[i] = 1;
rthreshold_converged = false; rthreshold_converged = false;
} }
rthreshold[i] = thisthreshold; rthreshold[i] = thisthreshold;
} }
} }
} }
com_x=0; com_x=0;
com_y=0; com_y=0;
//goto END; //goto END;
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++){
@ -361,21 +360,21 @@ 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;
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)
*/ */
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;
long np_counted = 0; long np_counted = 0;
float fbgr; float fbgr;
float backgroundMaxI=0; float backgroundMaxI=0;
@ -474,7 +473,6 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
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
*/ */
@ -538,7 +536,7 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
} }
} }
//END: ; //END: ;
free(temp); free(temp);
free(inx); free(inx);
@ -552,15 +550,14 @@ int peakfinder8(tPeakList *peaklist, float *data, char *mask, float *pix_r, long
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);
/*************************************************/ /*************************************************/
} }

4
peakfinder8/peakfinders.h
View File

@ -13,8 +13,8 @@
typedef struct { typedef struct {
public: public:
long nPeaks; long nPeaks;
long nHot; long nHot;
float peakResolution; // Radius of 80% of peaks float peakResolution; // Radius of 80% of peaks
float peakResolutionA; // Radius of 80% of peaks float peakResolutionA; // Radius of 80% of peaks
float peakDensity; // Density of peaks within this 80% figure float peakDensity; // Density of peaks within this 80% figure