From: Nate Rodning <rodning@relay.phys.ualberta.ca>
Date: Mon, 4 Oct 1999 09:38:30 -0600 (MDT)
To: e614tn@relay.phys.ualberta.ca
Subject: TN-38 - Laser Alignment System Studies
The attached technical note summarizes feasibility studies accomplished
thus far by George Price on the laser alignment system.
Please note that the appendices A1-A5 (containing C source code) are attached
as text files and appear below "inline". Appendices B containing additional
figures are in the .ps file.
nate
*************************************************
Nate Rodning, Professor of Physics
Associate Chair, Physics
University of Alberta
Edmonton, Alberta T6G 2N5
(780) 492 - 3518
Fax: 492 - 0714
Filename: alignmentsys.ps
Directory of Figures
/* qcam.c -- routines for accessing the Connectix QuickCam */
/* Version 0.4b, May 21, 1996 */
/******************************************************************
Copyright (C) 1996 by Scott Laird
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.
******************************************************************/
#include
#include
#include
#include
#include "qcam.h"
#include "qcamip.h"
/* filters() is a function added to the original source code to find the */
/* center of the diffraction patter captured by the camera */
void filters(scanbuf *pic8,struct qcam *q, float *cent, float *dataf);
/* vector() allocates space for (nh-nl) float data-values */
/* Filtered image data is read into the data location specified by vector()*/
float *vector(long nl, long nh);
/* moment() findes the average pixel intensity and standard deviaton of a picture. */
/* The moment() function is not used during normal camera operation. */
void moment(unsigned char *data, struct qcam *q);
void usage(void)
{
fprintf(stderr,"Usage:\n");
fprintf(stderr," qcam [options]\n");
fprintf(stderr," Options:\n");
fprintf(stderr," -x width Set width\n");
fprintf(stderr," -y height Set height\n");
fprintf(stderr," -p port Set port\n");
fprintf(stderr," -B bpp Set bits per pixel\n");
fprintf(stderr," -c val Set contrast\n");
fprintf(stderr," -w val Set white balance\n");
fprintf(stderr," -W Auto-set white balance\n");
fprintf(stderr," -b val Set brightness\n");
fprintf(stderr," -E \"vals\" Autoexposure mode, parameters required\n");
fprintf(stderr," -D Remove dark speckling\n");
fprintf(stderr," -H Display Histogram\n");
fprintf(stderr," -s val Set scaling factor (1, 2, or 4)\n");
fprintf(stderr," -t val Set top line of scan\n");
fprintf(stderr," -l val Set left column of scan\n");
fprintf(stderr," -f file Select configuration file\n");
fprintf(stderr," -u Force unidirectional mode\n");
fprintf(stderr," -V Show version information\n");
fprintf(stderr," -C Centroid calibrate mode\n");
fprintf(stderr," -F Size of low-pass filter\n");
}
void modegripe(struct qcam *q)
{
fprintf(stderr,"Unsupported resolution/bit depth!\n");
fprintf(stderr,"This program only supports 320x240, 160x120, and 80x60.\n");
fprintf(stderr,"You specified %d x %d. Try again.\n",q->width,q->height);
exit(1);
}
int main(int argc, char **argv)
{
scanbuf *scan = NULL;
int arg;
extern char *optarg;
extern int optind;
struct qcam *q;
int width, height;
int autocalibrate = 0;
int autoexposure = 0;
int ae_arg_cnt;
int ae_mode, ae_lum_target, ae_lum_tolerance;
int ae_lum_std_target, ae_lum_std_tolerance;
int dark = 0;
int histogram = 0;
int snapshots = 0;
int maxsnaps = 11;
static char *getargs = "hx:y:p:b:B:c:w:s:t:l:f:F:E:C:HVWuD";
int found;
float *centroid, *fcent;
long l=0,h=9,wth,i;
centroid=vector(l,h);
#ifdef DEBUG
struct timeval start, end;
double microseconds;
#endif
if(geteuid()) {
fprintf(stderr,"%s: Not installed SUID or run as root. Exiting.\n",
argv[0]);
exit(1);
}
q=qc_init();
/* unpleasant pre-search for -f option */
found = 0;
while ((arg = getopt(argc, argv, getargs)) > 0) {
if (arg == 'f' && !found) {
qc_initfile(q, optarg);
found = 1;
}
}
if (!found) qc_initfile(q, 0); /* use default */
qc_getresolution(q,&width,&height);
/* Read command line */
optind = 0;
while((arg=getopt(argc,argv,getargs))>0) {
switch (arg) {
case 'x':
width=atoi(optarg);
break;
case 'y':
height=atoi(optarg);
break;
case 'p':
if (!getuid())
q->port=strtol(optarg,NULL,0);
break;
case 'B':
qc_setbitdepth(q,atoi(optarg));
break;
case 'b':
qc_setbrightness(q,atoi(optarg));
break;
case 'c':
qc_setcontrast(q,atoi(optarg));
break;
case 'C':
qc_setcalibrate(q,atoi(optarg));
break;
case 'F':
qc_setfilter(q,atoi(optarg));
break;
case 'w':
qc_setwhitebal(q,atoi(optarg));
break;
case 's':
if (qc_settransfer_scale(q, atoi(optarg)))
fprintf(stderr, "Bad scaling factor, valid values are 1, 2 or 4\n");
break;
case 't':
if (qc_settop(q, atoi(optarg)))
fprintf(stderr, "Bad top line, valid values are 1 - 243\n");
break;
case 'l':
if (qc_setleft(q, atoi(optarg)))
fprintf(stderr, "Bad left column, valid values are 2 - 336, and must be even\n");
break;
case 'f':
/* already handled */
break;
case 'V':
fprintf(stderr,"%s: Version 0.8\n",argv[0]);
exit(0);
break;
case 'W':
autocalibrate = 1;
break;
case 'E':
autoexposure = 1;
ae_arg_cnt = sscanf(optarg, "%d %d %d %d %d", &ae_mode,
&ae_lum_target, &ae_lum_tolerance,
&ae_lum_std_target, &ae_lum_std_tolerance);
switch (ae_arg_cnt) {
case 5:
if (qcip_set_luminance_std_tolerance(q, ae_lum_std_tolerance))
fprintf(stderr, "Invalid luminance std tolerance.\n");
case 4:
if (qcip_set_luminance_std_target(q, ae_lum_std_target))
fprintf(stderr, "Invalid luminance std target.\n");
case 3:
if (qcip_set_luminance_tolerance(q, ae_lum_tolerance))
fprintf(stderr, "Invalid luminance tolerance.\n");
case 2:
if (qcip_set_luminance_target(q, ae_lum_target))
fprintf(stderr, "Invalid luminance target.\n");
case 1:
if (qcip_set_autoexposure_mode(ae_mode))
fprintf(stderr, "Invalid autoexposure mode.\n");
break;
default:
fprintf(stderr, "Invalid arguments for auto exposure.\n");
exit(1);
break;
}
break;
case 'H':
histogram = 1;
break;
case 'u':
qc_forceunidir(q);
break;
case 'D':
dark = 1;
break;
case 'h':
usage();
exit(0);
break;
default:
fprintf(stderr,"%s: Unknown option or error in option\n",argv[0]);
usage();
exit(1);
break;
}
}
if(qc_setresolution(q,width,height)) {
fprintf(stderr,"Invalid resolution. Exiting.\n");
exit(1);
}
/* Attempt to get permission to access IO ports. Must be root */
if (qc_open(q)) {
fprintf(stderr,"Cannot open QuickCam; exiting.\n");
exit(1);
}
setuid(getuid());
do {
/* Program the QuickCam */
qc_set(q);
if (autocalibrate) {
qc_calibrate(q);
}
/* Scan one image */
#ifdef DEBUG
if (gettimeofday(&start, NULL))
perror("get start time");
#endif
scan=qc_scan(q);
if (dark) {
fixdark(q, scan);
}
#ifdef DEBUG
if (gettimeofday(&end, NULL))
perror("get end time");
microseconds = end.tv_usec - start.tv_usec;
microseconds += 1000000.0 * (end.tv_sec - start.tv_sec);
fprintf(stderr, "%.0f ms to scan image", microseconds/1000.0);
if (microseconds != 0)
fprintf(stderr, " (%.1f fps)\n", 1000000.0/microseconds);
else
fprintf(stderr, "\n");
#endif
snapshots++;
} while (autoexposure && qcip_autoexposure(q, scan) && snapshots < maxsnaps);
/* Output the image to stdout */
if (histogram)
qcip_display_histogram(q, scan);
/* "wth" specifies the total number of pixel in filtered image */
/* "fcent" allocates memory for (wth x float) data values */
wth=(q->width*256)/(320*q->transfer_scale);
fcent=vector(0,wth*wth-1);
/* filters() finds centroid of intensity distribution */
/* x-centroid=centroid[0] */
/* y-centroid=centroid[1] */
/* error in x-centroid=centroid[2] */
/* error in y-centroid=centroid[3] */
/* sum of all pixel intensties for filered image = centroid[4] */
/* sum of all pixel intensities for non-filtered image = centroid[5] */
filters(scan,q,centroid,fcent);
/* Here the sum of all pixel intensities for the non-filtered image is found. */
for(i=0;i<=(q->width*q->height/(q->transfer_scale*q->transfer_scale))-1;i++) {
centroid[9]+=scan[i];
}
/* For testing and debugging all calculate values was piped to terminal */
printf("%f %f %f %f %f %f %f %f %f %f\n",centroid[0],centroid[1],centroid[2],centroid[3],centroid[4],centroid[5],centroid[6],centroid[7],centroid[8],centroid[9]);
/* The moment() function is used to find pixel intensity average and standard deviation */
/* Its disabled for normal operation. */
/* moment(scan,q);*/
/* writes non-filtered data to ouput file */
/* qc_writepgm(q,stdout,scan); */
free(scan);
/* Release IO privileges */
qc_close(q);
return 0;
}
#include
#include "math.h"
#include "include.h"
#include "filter.h"
#include "qcam.h"
/*Global Variables*/
long w,h;
long i, j, k;
unsigned char c;
FILE *out;
/* filters() finds the DFT or Discrete-Fourier-Transform of the raw image data */
/* then applies a low-pass filter to remove "high-frequency-noise". Finally, */
/* the centroid, error in centroid, and total-pixel-intensity-sum are found. */
void filters(unsigned char *pic8,struct qcam *q, float *cent, float *dataf)
{
float pi=3.14159265359;
double X1,X2,Xn1,Xn2,Y1,Y2,Yn1,Yn2,sin21,cos21,sin22,cos22,sincos1,sincos2,ex,ei;
float **spec,***array2;
float realx,imagx,realy,imagy,absy;
double phasex,phasey,errorx,errory;
int fltr;
/* The DFT is implimented with a FFT or Fast-Fourier-Transform algorithm */
/* FFT requires a window of dimension 2^N where "N" is an integer. The window */
/* size that the FFT is applied to is given by "w" and "h". All pixel outside */
/* of the w by h window are set to zero. */
/* Choose the size of window to filter*/
w=(q->width*256)/(320*q->transfer_scale);
h=w;
/* Make 2-d array from 1-d array produced by qcam. I needed to do this */
/* the code used to find the FFT would work properly. */
array2=f3tensor(1,1,1,h,1,w);
/*read data into 2-d array*/
for(i=1;i<=h;i++) {
for(j=1+w/8;j<=w+w/8;j++) {
if(ih-h/32-1) array2[1][i][j-w/8]= 0.0;
else array2[1][i][j-w/8]=pic8[(i-h/32-1)*(q->width/q->transfer_scale)+j-1];
}
}
/*initialize matrix for use with rflt3(). rflt3() is the function which */
/*finds the DFT using FFT. This function is in a seperate file labeled rflt3.c. */
/*This matrix provides space for the N/2+1 freqency componets of DFT. However, I've */
/*made all high freqency componets equal to zero so the matrix below is set to zero. */
/*But the rflt3() still wants this matrix (I think its useful when considering aliasing). */
spec=matrix(1,1,1,2*h);
/*find DFT using rlft3()*/
rlft3(array2,spec,1,h,w,1);
/*find X-section of DFT */
/* for(i=1;i<=h;i++){
absy=sqrt(array2[1][i][1]*array2[1][i][1]+array2[1][i][2]*array2[1][i][2]);
printf("%f ",absy);
} */
/*find phase of array[1][1][1]*/
/*phasex is phase of x-direction-harmonic converted to pixels*/
/*phasey is phase of y-direction-harmonic converted to pixels*/
realx=array2[1][1][3];
imagx=array2[1][1][4];
realy=array2[1][2][1];
imagy=array2[1][2][2];
/* for(i=1;i<=h;i++){
for(j=1;j<=w;j++){
realx+=array2[1][i][j]*cos(2*pi*j/w);
imagx+=array2[1][i][j]*sin(2*pi*j/w);
realy+=array2[1][i][j]*cos(2*pi*i/h);
imagy+=array2[1][i][j]*sin(2*pi*i/h);
}
} */
phasex=atan(imagx/realx);
phasex=h*phasex/6.28318530717959;
phasey=atan(imagy/realy);
phasey=h*phasey/6.28318530717959;
cent[0]=phasex;
cent[1]=phasey;
/*find error in x-centroid and y-centroid using phase shift method*/
/*errorx=error in phase for x-direction-harmonic converted to pixels.*/
/*errory=error in phase for y-direction-harmonic converted to pixels.*/
/*ex=error in x-position of pixel*/
/*ey=error in y-postiion of pixel*/
for(i=1;i<=h;i++) {
sin21+=w*sin(2*pi*j/w)*sin(2*pi*j/w);
cos21+=w*cos(2*pi*j/w)*cos(2*pi*j/w);
sin22+=h*sin(2*pi*i/h)*sin(2*pi*i/h);
cos22+=h*cos(2*pi*i/h)*cos(2*pi*i/h);
}
ex=0.02;
ei=1.0;
/*errorx=sqrt((w*w*ei*ei*)/((4*pi*pi)(X1*X1+Y1*Y1)) + ex*ex);*/
/*errory=sqrt(((h*h*ei*ei)/((4*pi*pi)(X2*X2+Y2*Y2))) + ex*ex);*/
errorx=sqrt(((w*w*ei*ei)*(sin21+cos21)/((4*pi*pi)*(realx*realx+imagx*imagx))) + ex*ex);
errory=sqrt(((w*w*ei*ei)*(sin22+cos22)/((4*pi*pi)*(realy*realy+imagy*imagy))) + ex*ex);
cent[2]=errorx;
cent[3]=errory;
/* printf("%1f %1f errorx=%1f errory=%1f X1=%f Y1=%f sin21=%f cos21=%f realx=%f imagx=%f", phasex, phasey, errorx, errory, X1, Y1, sin21, cos21, realx, imagx);*/
/*printf("%1f %1f ",phasex,phasey);*/
if(q->cal==1){
/* Apply low-pass filter. */
/* This function is located in Filter-lib.c */
fltr=q->filter;
lowpass(array2,spec,fltr,w,h);
/*find Inverse Discrete Fourier Transform or IDFT of filtered data*/
rlft3(array2,spec,1,h,w,-1);
/*write filtered data to 1-d array for use by TCL */
for(i=1;i<=h;i++) {
for(j=1;j<=w;j++) {
dataf[(i-1)*h+(j-1)]=array2[1][i][j];
}
}
/* write array2 data to output file */
if((out = fopen("dataout.pgm","w")) ==NULL){
printf("error opening file\n");
exit(1);
}
fprintf(out,"P5\n");
fprintf(out,"%d %d\n",w,h);
fprintf(out,"%d\n",63);
for(i=0;i<=w*h-1;i++) {
c=dataf[i];
fputc(c,out);
}
close(out);
/* find the centroid of filtered data. fcentroid() is a function */
/* located in a file "fcentroid.c" */
fcentroid(array2,q,w,h,cent);
}
}
#include
#include
#include
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[]);
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl...nrh][ncl...nch] */
{
long i, nrow=nrh-nrl+1, ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with rnage t[nrl...nrh][ncl...nch][ndl...ndh] */
{
long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1, ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
void lowpass(float ***array2, float ** speq, int filter, long w, long h)
/* apply lowpass filter to discrete fourier transfor of image*/
{
long i, j, k;
/* make FFT of image "zero" except for difined freqency window */
for(i=1;i filter*filter/4
&& (h-i)*(h-i)+j*j/4 > filter*filter/4) array2[1][i][j]=0.0;
else array2[1][i][j]=(2.0/(w*h))*array2[1][i][j];
}
}
/*make speq=0*/
for(i=1;i<2*h+1;i++) {
speq[1][i]=0.0;
}
}
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n", error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
/* Allocate a float vector with subscript range v[nl..nh]*/
float *vector(long nl, long nh)
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
return v-nl+NR_END;
}
#include
void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2,
unsigned long nn3, int isign)
{
void fourn(float data[], unsigned long nn[], int ndim, int isign);
void nrerror(char error_text[]);
unsigned long i1,i2,i3,j1,j2,j3,nn[4],ii3;
double theta,wi,wpi,wpr,wr,wtemp;
float c1,c2,h1r,h1i,h2r,h2i;
if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3)
nrerror("rlft3: problem with dimensions or contiguity of data array\n");
c1=0.5;
c2 = -0.5*isign;
theta=isign*(6.28318530717959/nn3);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
nn[1]=nn1;
nn[2]=nn2;
nn[3]=nn3 >> 1;
if (isign == 1) {
fourn(&data[1][1][1]-1,nn,3,isign);
for (i1=1;i1<=nn1;i1++)
for (i2=1,j2=0;i2<=nn2;i2++) {
speq[i1][++j2]=data[i1][i2][1];
speq[i1][++j2]=data[i1][i2][2];
}
}
for (i1=1;i1<=nn1;i1++) {
j1=(i1 != 1 ? nn1-i1+2 : 1);
wr=1.0;
wi=0.0;
for (ii3=1,i3=1;i3<=(nn3>>2)+1;i3++,ii3+=2) {
for (i2=1;i2<=nn2;i2++) {
if (i3 == 1) {
j2=(i2 != 1 ? ((nn2-i2)<<1)+3 : 1);
h1r=c1*(data[i1][i2][1]+speq[j1][j2]);
h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]);
h2i=c2*(data[i1][i2][1]-speq[j1][j2]);
h2r= -c2*(data[i1][i2][2]+speq[j1][j2+1]);
data[i1][i2][1]=h1r+h2r;
data[i1][i2][2]=h1i+h2i;
speq[j1][j2]=h1r-h2r;
speq[j1][j2+1]=h2i-h1i;
} else {
j2=(i2 != 1 ? nn2-i2+2 : 1);
j3=nn3+3-(i3<<1);
h1r=c1*(data[i1][i2][ii3]+data[j1][j2][j3]);
h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]);
h2i=c2*(data[i1][i2][ii3]-data[j1][j2][j3]);
h2r= -c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]);
data[i1][i2][ii3]=h1r+wr*h2r-wi*h2i;
data[i1][i2][ii3+1]=h1i+wr*h2i+wi*h2r;
data[j1][j2][j3]=h1r-wr*h2r+wi*h2i;
data[j1][j2][j3+1]= -h1i+wr*h2i+wi*h2r;
}
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
}
if (isign == -1)
fourn(&data[1][1][1]-1,nn,3,isign);
}
#include
#include "math.h"
#include "include.h"
#include "qcam.h"
/*fcentroid() finds the center and error in center of the diffraction pattern located in "data"*/
/*All calculated data is stored in pointer-array "cen"*/
void fcentroid(float ***data, struct qcam *q, long w, long h, float *cen)
{
long dy, dx, i, j, x, y, xmax, xmin, ymax, ymin ,n ,count;
float c, totcnt, row,col, maxr,maxc, min, total1, total ,total2 , row1, col1, out[5];
float xcen,ycen;
float yerr, xerr, yr, xr;
float sumx,sumy,sumI,sumI2,sumx2,sumy2,sumxI,sumyI;
/* find maximum counts in rows and columns */
col1=0.0;
row1=0.0;
/* find row with greatest pixel sum */
for(i=1;i<=h;i++) {
for(j=1;j<=w;j++) {
row+=data[1][i][j];
}
if(row>row1) {
maxr=row;
y=i;
}
total1+=row;
row1=maxr;
row=0.0;
}
cen[8]=total1;
/* find column with greatest pixel sum */
for(j=1;j<=w;j++) {
for(i=1;i<=h;i++) {
col+=data[1][i][j];
}
if(col>col1) {
maxc=col;
x=j;
}
total2+=col;
col1=maxc;
col=0.0;
}
/* find size of pixel-window to find centroid of */
if(y<=h/2) {
ymax=2*y-1;
ymin=1;
}
else {
ymax=h;
ymin=2*y-h;
}
if(x<=w/2) {
xmax=2*x-1;
xmin=1;
}
else {
xmax=w;
xmin=2*x-w;
}
dx=xmax-xmin;
dy=ymax-ymin;
if(dy<=dx) dx=dy;
else dy=dx;
/* find x and y coordinated of diffraction pattern center */
/* "xcen"= X-coordintat of diffraction pattern center */
/* "ycen"= Y-coordinate of diffraction pattern center */
for(i=y-dy/2;i<=y+dy/2;i++) {
for(j=x-dx/2;j<=x+dx/2;j++) {
xcen+=j*data[1][i][j];
ycen+=i*data[1][i][j];
total+=data[1][i][j];
/* find error in "x" and "y" coordinated of diffraction pattern center. */
/* "xerr"= Error in X-position */
/* "yerr"= Error in Y-position */
sumI+=data[1][i][j];
sumI2+=data[1][i][j]*data[1][i][j];
sumx+=j;
sumy+=i;
sumx2+=j*j;
sumy2+=i*i;
sumxI+=j*data[1][i][j];
sumyI+=i*data[1][i][j];
}
}
xcen=xcen/total;
cen[4]=xcen;
ycen=ycen/total;
cen[5]=ycen;
xr=0.5;
yr=1.0;
n=(xmax-xmin+1)*(ymax-ymin+1);
xerr=sqrt((xr*xr)*sumI2/(sumI*sumI)+(yr*yr)*((sumI*sumI)*(sumx2)-2*(sumxI*sumI*sumx)+n*(sumxI*sumxI))/(sumI*sumI*sumI*sumI));
yerr=sqrt((xr*xr)*sumI2/(sumI*sumI)+(yr*yr)*((sumI*sumI)*(sumy2)-2*(sumyI*sumI*sumy)+n*(sumyI*sumyI))/(sumI*sumI*sumI*sumI));
if(sumI<1) {
xerr=0.0;
yerr=0.0;
}
cen[6]=xerr;
cen[7]=yerr;
/* printf("x=%d, y=%d, xcen=%f, ycen=%f, xerr=%f, yerr=%f, total counts=%f\n",x ,y ,xcen,ycen,xerr,yerr,total);*/
/* printf("sumI2=%f, sumI=%f, xmax=%d, xmin=%d, ymax=%d, ymin=%d\n",sumI2,sumI,xmax,xmin,ymax,ymin);*/
/* find histogram */
/* for(c=0;c<=63;c++) {
count=0;
for(i=ymin;i<=ymax;i++) {
for(j=xmin;j<=xmax;j++) {
if(c==data[1][i][j]) count=count++;
}
}
printf("c=%f count=%ld \n",c,count);
}*/
}
TN-38 / Nate Rodning
- Created for the The Center for Subatomic Research E614 Project Projects Page.
- Created by The CoCoBoard.