pshell_config/dev
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Startup

Browse Source
This commit is contained in:
gobbo_a
2019-08-16 09:54:45 +02:00
parent c24eb95fcd
commit 444cfe9c6e
4 changed files with 636 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
script/cpython/image_functions.py Normal file
View File

@@ -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)

550
script/ms/image_get_int.c Normal file
View File

@@ -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 <filter_nsigma>: 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 <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <time.h>
/* 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 <image_info_filename> <thresh1> <thresh2> <thresh3> <thresh4> [<-median> <filter_nsigma>]\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<PIXEL_TOTAL; i++);
{
img[i] = (unsigned long int) img16[i];
}
}
else if(image_color_depth == 32)
{
pixel_count = fread(img,sizeof(unsigned long int),
PIXEL_TOTAL,file_img);
}
else
{
fprintf(stderr,"Color-depth %d not supported.\n",image_color_depth);
return 1;
}
fclose(file_img);
retry++;
wait(0.01);
}
if (pixel_count != PIXEL_TOTAL)
{
fprintf(stderr,"Can't read pixel data.\n");
fprintf(stderr,"read %d pixels of %d\n",pixel_count,PIXEL_TOTAL);
return 1;
}
/*-------------------------------------*/
/* optionally perform median filtering */
if (filter_nsigma > 0)
{
unsigned long int filt_img[PIXEL_TOTAL];
/* filter only bounding box of ROIs */
int fx1 = (bx1<x1) ? bx1 : x1;
int fy1 = (by1<y1) ? by1 : y1;
int fx2 = (bx2>x2) ? 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 = (bx1<x1) ? bx1 : x1;
int fy1 = (by1<y1) ? by1 : y1;
int fx2 = (bx2>x2) ? 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 $ ===*/

57
script/ms/image_get_int.py Normal file
View File

@@ -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))

3
script/ms/test.py Normal file
View File

@@ -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)