diff --git a/script/cpython/image_functions.py b/script/cpython/image_functions.py new file mode 100644 index 0000000..dab5206 --- /dev/null +++ b/script/cpython/image_functions.py @@ -0,0 +1,26 @@ +import numpy + +def img_get_int(fname, thres1, thres2, thres3, thres4, header, width, height, depth, x1,y1,x2,y2, bx1,by1,bx2,by2 , filter_median = False, filter_nsigma = 0): + # read actual image file + img = numpy.fromfile(fname, dtype=numpy.uint32) + img.shape = height, width + # signal roi + area_I = ( x2 - x1 + 1) * ( y2 - y1 + 1) + I_sum = img[y1:y2, x1:x2].sum() + thresh1_count = len(numpy.where(img>thres1)[0]) + thresh2_count = len(numpy.where(img>thres2)[0]) + thresh3_count = len(numpy.where(img>thres3)[0]) + thresh4_count = len(numpy.where(img>thres4)[0]) + + # background roi + I_sum_bgr = img[by1:by2, bx1:bx2].sum() + area_bgr= (bx2 - bx1 + 1) * (by2 - by1 + 1) + return (I_sum, area_I, thresh1_count, thresh2_count, thresh3_count, thresh4_count, I_sum_bgr, area_bgr) + + +def img_read(fname, header, width, height, depth): + img = numpy.fromfile(fname, dtype=numpy.uint32) + img.shape = height, width + #return img.flatten().astype(int) + return img.astype(int) + \ No newline at end of file diff --git a/script/ms/image_get_int.c b/script/ms/image_get_int.c new file mode 100644 index 0000000..927f291 --- /dev/null +++ b/script/ms/image_get_int.c @@ -0,0 +1,550 @@ +/*============================================================================= + +Program name: image_get_int.c + =============== + +$Date: 2009/05/27 15:13:48 $ +$Author: schlepuetz $ +$Source: /cvs/G/SPEC/local/X04SA/ES3/src/image_get_int.c,v $ +$Revision: 1.4 $ +$Tag: $ + +Description: This program extracts various intensity information from binary + image files. + +Note: This macro is based on an earlier version by O. Bunk, + November 2004. + +Authors(s): O. Bunk (OB), C. M. Schlepuetz (CS) +Co-author(s): R. Herger (RH), P. R. Willmott (PW) +Address: Surface Diffraction Station + Materials Science Beamline X04SA + Swiss Light Source (SLS) + Paul Scherrer Institute + CH - 5232 Villigen PSI +Created: 2005/07/13 + +Compile: g++ -o image_get_int image_get_int.c + +Change Log: +----------- + +2005/07/13 (CS): +- added information and comments, changed formatting to comply with X04SA + conventions. + +2005/12/12 (CS): +- implemented flexibility to read different binary file formats (different + image sizes, header lengths, color depths). +- implemented a fourth threshold level. +- information on the image which is to be evaluated (name, header length, + dimensions, box coordinates, etc.) is now stored in a file which is passed as + a command line argument. + +2006/01/27 (CS): +- fixed wrong order of arguments for median3x3 + +2006/09/21 (CS,PW): +- created mode for selectively median filtering of random hot and dead pixels: + - additional input argument : specifies by how many standard + deviations the value of a pixel is allowed to differ from the average of + its surrounding pixels. + +2006/11/29 (PW,CS): +- added the wait() function. +- included wait loop when reading the image file to check that the file is + complete. If file is not complete after 10 seconds (1000 * 0.01 seconds), an + error message is issued and the program aborted. + +=============================================================================*/ + +#include +#include +#include +#include +#include + +/* uncomment the next lines to print debug info during execution */ +//#define print_debug_info +//#define print_debug_info_2 +//#define print_debug_info_3 + +#define PIXEL_TOTAL (image_xdim * image_ydim) +#define I(x,y) (img[(x) + image_xdim*(y)]) +#define FALSE 0 + +/*---------------------------------------------------------------------------*/ +/* wait function */ +/*---------------------------------------------------------------------------*/ + +void wait ( double seconds ) +{ + clock_t endwait; + endwait = clock () + int(floor(seconds * CLOCKS_PER_SEC)) ; + while (clock() < endwait) {} +} + +/*---------------------------------------------------------------------------*/ +/* functions to perform a 3x3 median filtering */ +/*---------------------------------------------------------------------------*/ + +#define MED_BOX 3 +#define MED_BOX_HALF (MED_BOX/2) + +/*---------------------------------------------------------------------------*/ +/* sub-function for the median filter below */ +/* needed to sort the array elements */ +/* refer to the qsort documentation in cstdlib.h */ + +int qsort_comp_func(const void *I1,const void *I2) +{ + if (*(int *)I1 < *(int *)I2) + return -1; + if (*(int* )I1 > *(int *)I2) + return 1; + return 0; +} + +/*---------------------------------------------------------------------------*/ +/* 3x3 median filter to the image data *img */ + +void median_3x3(unsigned long int *img,unsigned long int *filt_img,int image_xdim,int image_ydim,int x1,int y1,int x2,int y2) +{ + int x,y; + + if (x1 < MED_BOX_HALF) + x1 = MED_BOX_HALF; + if (y1 < MED_BOX_HALF) + y1 = MED_BOX_HALF; + if (x2 > image_xdim -1 -MED_BOX_HALF) + x2 = image_xdim -1 -MED_BOX_HALF; + if (y2 > image_ydim -1 -MED_BOX_HALF) + y2 = image_ydim -1 -MED_BOX_HALF; + #ifdef print_debug_info + printf("median: (x1,y1)=(%d %d) (x2,y2)=(%d %d) box-size=%d half-box %d\n", + x1,y1,x2,y2,MED_BOX,MED_BOX_HALF); + #endif + + for (x=x1;x <= x2;x++) + for (y=y1;y <= y2;y++) + { + int I_values[MED_BOX*MED_BOX]; + int xm,ym,i = 0; + for (xm=x-MED_BOX_HALF;xm <= x+MED_BOX_HALF;xm++) + for (ym=y-MED_BOX_HALF;ym <= y+MED_BOX_HALF;ym++) + { + int img_index = xm + ym*image_xdim; + +#ifdef print_debug_info + printf("median box %5d %5d %5d\n",xm,ym,img[img_index]); +#endif + + I_values[i++] = img[img_index]; + } + qsort(I_values,i,sizeof(*I_values),qsort_comp_func); + +#ifdef print_debug_info + for (i=0;i < MED_BOX*MED_BOX;i++) + printf("qsort %d %5d\n",i,I_values[i]); + printf("%d %d: %d\n",x,y,I_values[(MED_BOX*MED_BOX)/2]); +#endif + filt_img[x + y*image_xdim] = I_values[(MED_BOX*MED_BOX)/2]; + } + + for (x=x1;x <= x2;x++) + for (y=y1;y <= y2;y++) + { + int img_index = x + y*image_xdim; + img[img_index] = filt_img[img_index]; + } + return; +} + +/*---------------------------------------------------------------------------*/ +/* main function image_get_int */ +/* ============= */ +/*---------------------------------------------------------------------------*/ + +int main(const int no_of_params,const char **params) +{ + FILE *file_img_info; + char filename_img[255]; + FILE *file_img; + double image_timestamp; + int image_header_length; + int image_xdim,image_ydim,image_color_depth; + int x1,y1,x2,y2; + int bx1,by1,bx2,by2; + int x,y,i; + long int I_sum = 0,I_sum_bgr = 0,area_I,area_bgr; + int thresh1_count = 0,thresh2_count = 0,thresh3_count = 0,thresh4_count = 0; + int filter_median = 0; + int parameter_okay = 1; + int thresh1,thresh2,thresh3,thresh4; + int retry; + float filter_nsigma; + + /* check input parameter */ + if ((no_of_params < 6) || + (no_of_params > 8) || + (sscanf(params[ 2],"%d",&thresh1) != 1) || + (sscanf(params[ 3],"%d",&thresh2) != 1) || + (sscanf(params[ 4],"%d",&thresh3) != 1) || + (sscanf(params[ 5],"%d",&thresh4) != 1) || + (thresh1 < 0) || (thresh2 < 0) || (thresh3 < 0) || (thresh4 < 0)){ + parameter_okay = FALSE; + } + +#ifdef print_debug_info + printf("No. of command-line parameters: %d\n",no_of_params); +#endif + + if (no_of_params > 6) + { + if (strcmp(params[ 6],"-median") == 0) + filter_median = !FALSE; + else + parameter_okay = FALSE; + } + + if (no_of_params > 7) + { + +#ifdef print_debug_info_2 + printf("params(7) = %s\n",params[7]); +#endif + + if (sscanf(params[ 7],"%g",&filter_nsigma) != 1) + { + parameter_okay = FALSE; + } + if (filter_nsigma < 0) + { + parameter_okay = FALSE; + } + +#ifdef print_debug_info_2 + printf("filter_nsigma = %f\n",filter_nsigma); +#endif + + } + else + { + filter_nsigma = 0; + } + + /* open image info text-file for read-only access */ + if ((file_img_info = fopen(params[1],"r")) == NULL) + { + fprintf(stderr,"Can't open image info file %s\n",params[1]); + return 1; + } + else + { + +#ifdef print_debug_info + printf("read image info file %s\n",params[1]); +#endif + + int read_count = fscanf (file_img_info, + "%s %f %d %d %d %d %d %d %d %d %d %d %d %d", + filename_img,&image_timestamp,&image_header_length, + &image_xdim,&image_ydim,&image_color_depth, + &x1,&y1,&x2,&y2, + &bx1,&by1,&bx2,&by2); + + if (read_count != 14) + { + fprintf(stderr,"Can't read image info file %s\n",params[1]); + return 1; + } + fclose(file_img_info); + + /* check consistency of image information */ + if ((x1 < 0) || (y1 < 0) || (x2 < 0) || (y2 < 0) || + (x1 >= image_xdim) || (y1 >= image_ydim) || + (x2 >= image_xdim) || (y2 >= image_ydim) || + (x2 < x1) || (y2 < y1) || + (bx1 < 0) || (by1 < 0) || (bx2 < 0) || (by2 < 0) || + (bx1 >= image_xdim) || (by1 >= image_ydim) || + (bx2 >= image_xdim) || (by2 >= image_ydim) || + (bx2 < bx1) || (by2 < by1)) + { + parameter_okay = FALSE; + } + } + + /* print help text in case of wrong parameters */ + if (!parameter_okay) + { + fprintf(stderr,"Get the number of pixel within a box (x1,y1,x2,y2)\n"); + fprintf(stderr,"with counts above four threshold values.\n"); + fprintf(stderr,"(0,0) is the upper left corner of the image.\n"); + fprintf(stderr,"(%d,%d) is the lower right corner.\n", + image_xdim-1,image_ydim-1); + fprintf(stderr,"Median filtering can be used to remove dead and hot pixels.\n\n"); + fprintf(stderr,"A second integration box can be used to determine the background\n"); + fprintf(stderr,"(coordinates (bx1,by1) (bx2,by2)).\n"); + fprintf(stderr,"Usage:\n"); + fprintf(stderr,"%s [<-median> ]\n\n", + params[0]); + fprintf(stderr,"All parameters are mandatory except the last -median switch and filter_nsigma value.\n"); + fprintf(stderr,"If the filter_nsigma value is set to zero or is unspecified, all the image will be automatically median filtered.\n"); + fprintf(stderr,"Positive values of filter_nsigma will lead to selective median filtering of individual pixels that differ by more than (filter_nsigma x st. dev.) of the eight neighbouring pixels.\n"); + fprintf(stderr,"The order of the parameters is fixed.\n"); + fprintf(stderr,"The image information file is automatically created by SPEC.\n"); + return 1; + } + + unsigned long int img[PIXEL_TOTAL]; + + area_I = (x2- x1+1)*( y2- y1+1); + area_bgr = (bx2-bx1+1)*(by2-by1+1); + +#ifdef print_debug_info + printf("( x1, y1)=(%3d %3d) ( x2, y2)=(%3d %3d) area=%ld\n", + x1,y1,x2,y2,area_I); + printf("(bx1,by1)=(%3d %3d) (bx2,by2)=(%3d %3d) area=%ld\n", + bx1,by1,bx2,by2,area_bgr); + if (filter_median) + printf("with median filtering\n"); + else + printf("no median filtering\n"); +#endif + + /* open the actual image as text-file for read-only access */ + retry = 0; + int pixel_count = 0; + while ((pixel_count != PIXEL_TOTAL) & (retry < 1000)){ + + /*fprintf(stderr,"Waiting for file %s\tretry = %d\n",filename_img,retry);*/ + + if ((file_img = fopen(filename_img,"r")) == NULL) + { + fprintf(stderr,"Can't open %s\n",filename_img); + return 1; + } + + if(image_header_length >0) + { + fseek (file_img,image_header_length,SEEK_SET); + } + + if(image_color_depth == 16) + { + unsigned short int img16[PIXEL_TOTAL]; + pixel_count = fread(img16,sizeof(unsigned short int), + PIXEL_TOTAL,file_img); + for (i=0; i 0) + { + unsigned long int filt_img[PIXEL_TOTAL]; + + /* filter only bounding box of ROIs */ + int fx1 = (bx1x2) ? bx2 : x2; + int fy2 = (by2>y2) ? by2 : y2; + + /* make sure we do not operate on the image border */ + fx1 = (fx1>0) ? fx1 : 1; + fy1 = (fy1>0) ? fy1 : 1; + fx2 = (fx2<(image_xdim-1)) ? fx2 : (image_xdim-2); + fy2 = (fy2<(image_ydim-1)) ? fy2 : (image_ydim-2); + + for (x=fx1; x<=fx2; x++){ + for (y=fy1; y<=fy2; y++){ + + /* calculate average of eight surrounding pixels */ + int xb,yb,i = 0; + double box_mean = 0; + for (xb=x-MED_BOX_HALF;xb <= x+MED_BOX_HALF;xb++){ + for (yb=y-MED_BOX_HALF;yb <= y+MED_BOX_HALF;yb++){ + int img_index = xb + yb*image_xdim; + +#ifdef print_debug_info + printf("median box %5d %5d %5d\n",xb,yb,img[img_index]); +#endif + if (!((xb == x) && (yb == y))) + { + box_mean += img[img_index]; + } + } + } + box_mean = box_mean / (MED_BOX*MED_BOX-1); + +#ifdef print_debug_info_2 + printf("box_mean = %f\n",box_mean); +#endif + + /* calculate standard deviation of eight surrounding pixels */ + xb = 0; yb = 0; i = 0; + double box_sigma = 0; + for (xb=x-MED_BOX_HALF;xb <= x+MED_BOX_HALF;xb++){ + for (yb=y-MED_BOX_HALF;yb <= y+MED_BOX_HALF;yb++){ + int img_index = xb + yb*image_xdim; + +#ifdef print_debug_info + printf("median box %5d %5d %5d\n",xb,yb,img[img_index]); +#endif + if (!((xb == x) && (yb == y))) + { + box_sigma += pow((img[img_index]-box_mean),2); + } + } + } + box_sigma = sqrt(box_sigma / (MED_BOX*MED_BOX-1)); + +#ifdef print_debug_info_2 + printf("box_sigma = %f\n",box_sigma); + printf("filter_nsigma = %f\n",filter_nsigma); +#endif + + if (fabs(I(x,y)-box_mean) > (filter_nsigma * box_sigma)) + { + +#ifdef print_debug_info_3 + printf("filtered (%d,%d): %d --> ",x,y,I(x,y)); +#endif + + median_3x3(img,filt_img,image_xdim,image_ydim,x,y,x,y); + +#ifdef print_debug_info_3 + printf("%d\n",I(x,y)); +#endif + + } + } + } + } + else if (filter_median) + { + + unsigned long int filt_img[PIXEL_TOTAL]; + + /* filter only bounding box of ROIs */ + int fx1 = (bx1x2) ? bx2 : x2; + int fy2 = (by2>y2) ? by2 : y2; + + /* make sure we do not operate on the image border */ + fx1 = (fx1>0) ? fx1 : 1; + fy1 = (fy1>0) ? fy1 : 1; + fx2 = (fx2<(image_xdim-1)) ? fx2 : (image_xdim-2); + fy2 = (fy2<(image_ydim-1)) ? fy2 : (image_ydim-2); + + + /* perform median filtering on the whole image */ + median_3x3(img,filt_img,image_xdim,image_ydim,fx1,fy1,fx2,fy2); + } + + /* calculate the results (1 of 2): intensity integration area */ + for (x=x1;x <= x2;x++) + { + for (y=y1;y <= y2;y++) + { +#ifdef print_debug_info + printf("(%5d,%5d) / index %5d: %5u\n", + x,y,x+y*image_xdim,img[x+y*image_xdim]); +#endif + I_sum += I(x,y); + if (I(x,y) > thresh1) + thresh1_count++; + if (I(x,y) > thresh2) + thresh2_count++; + if (I(x,y) > thresh3) + thresh3_count++; + if (I(x,y) > thresh4) + thresh4_count++; + } + } + + /* calculate the results (2 of 2): background integration area */ + for (x=bx1;x <= bx2;x++) + { + for (y=by1;y <= by2;y++) + { + /* only sum up background outside the intensity integration box */ + if ((x < x1) || (x > x2) || (y < y1) || (y > y2)) + { +#ifdef print_debug_info + printf("(%5d,%5d) / index %5d: %5u\n", + x,y,x+y*image_xdim,img[x+y*image_xdim]); +#endif + I_sum_bgr += I(x,y); + } + /* if pixel is inside the intensity integration box subtract it */ + /* from the background box area */ + else + { + area_bgr += -1; + } + } + } + + /* print results */ + printf("I_sum %ld area_I %ld Threshold1 %d Threshold2 %d Threshold3 %d Threshold4 %d I_bgr %ld area_bgr %ld\n", + I_sum,area_I,thresh1_count,thresh2_count,thresh3_count,thresh4_count, + I_sum_bgr,area_bgr); + return 0; +} + +/*============================================================================ + +#################################################### + emacs setup: force text mode to get no help with # + indentation and force use of spaces # + when tabbing. # + Local Variables: # + mode:text # + indent-tabs-mode:nil # + End: # +#################################################### + +$Log: image_get_int.c,v $ +Revision 1.4 2009/05/27 15:13:48 schlepuetz +commented out status message when waiting for file on server + +Revision 1.3 2006/12/01 08:31:00 schlepuetz +added the wait() function; included wait loop when reading the image file to check that the file is complete + +Revision 1.2 2006/09/21 11:42:39 schlepuetz +added mode for selectively filtering random hot and dead pixels + +Revision 1.1 2006/03/02 13:24:27 schlepuetz +Program to extract intensities from pixel detector data - first release, verified to work on PILATUS I and PILATUS II in conjunction with .img files. + + +============================= End of $RCSfile: image_get_int.c,v $ ===*/ diff --git a/script/ms/image_get_int.py b/script/ms/image_get_int.py new file mode 100644 index 0000000..ebc8884 --- /dev/null +++ b/script/ms/image_get_int.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python +import sys + +# parse aguments +# image info file +# threshold 1 +# threshold 2 +# threshold 3 +# threshold 4 +# median filter +# filter nsigma +info = sys.argv[1] +thres1 = int(sys.argv[2]) +thres2 = int(sys.argv[3]) +thres3 = int(sys.argv[4]) +thres4 = int(sys.argv[5]) + +filter_median = False +if len(sys.argv) > 6 and sys.argv[6] == '-median': + filter_median = True + +filter_nsigma = 0.0 +if len(sys.argv) > 7: + filter_nsigma = float(sys.argv[7]) + +# read image info text file +f = open(info) + +fname = f.next().strip() +timestamp = f.next().strip() +header = f.next().strip() +width, height, depth = [int(p) for p in f.next().strip().split()] +x1,y1,x2,y2 = [int(p) for p in f.next().strip().split()] +bx1,by1,bx2,by2 = [int(p) for p in f.next().strip().split()] + +f.close() + +# read actual image file +import numpy +img = numpy.fromfile(fname, dtype=numpy.uint32) +img.shape = height, width + +# signal roi +area_I = ( x2 - x1 + 1) * ( y2 - y1 + 1) + +I_sum = img[y1:y2, x1:x2].sum() +thresh1_count = len(numpy.where(img>thres1)[0]) +thresh2_count = len(numpy.where(img>thres2)[0]) +thresh3_count = len(numpy.where(img>thres3)[0]) +thresh4_count = len(numpy.where(img>thres4)[0]) + +# background roi +I_sum_bgr = img[by1:by2, bx1:bx2].sum() +area_bgr= (bx2 - bx1 + 1) * (by2 - by1 + 1) + +set_return ((I_sum, area_I, thresh1_count, thresh2_count, thresh3_count, thresh4_count, I_sum_bgr, area_bgr)) + diff --git a/script/ms/test.py b/script/ms/test.py new file mode 100644 index 0000000..acdcfc3 --- /dev/null +++ b/script/ms/test.py @@ -0,0 +1,3 @@ +ret = img_get_int(tmp_file, threshold1, threshold2, threshold3, threshold4, \ + self.pixel.image_header_length, self.pixel.PIX_XDIM, self.pixel.PIX_YDIM,self.pixel.PIX_COLOR_DEPTH, \ + roi[0], roi[1], roi[2], roi[3], bkroi[0], bkroi[1], bkroi[2], bkroi[3], filter_median = True, filter_nsigma = 30)