551 lines
17 KiB
C
551 lines
17 KiB
C
/*=============================================================================
|
|
|
|
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 $ ===*/
|