aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/calibrate_detector.c213
1 files changed, 109 insertions, 104 deletions
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 25d7a3c8..54472b04 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -51,6 +51,103 @@ static void show_help(const char *s)
"\n");
}
+//FIXME: should this function be in detector.h?
+int calculate_projected_peak(struct panel *panel, struct rvec q,
+ double kk, double *fs, double *ss)
+{
+
+ if (panel == NULL) return 1;
+
+ double xd, yd, cl;
+ double plx, ply;
+
+ const double den = kk + q.w;
+
+ /* Camera length for this panel */
+ cl = panel->clen;
+
+ /* Coordinates of peak relative to central beam, in m */
+ xd = cl * q.u / den;
+ yd = cl * q.v / den;
+
+ /* Convert to pixels */
+ xd *= panel->res;
+ yd *= panel->res;
+
+ /* Convert to relative to the panel corner */
+ plx = xd - panel->cnx;
+ ply = yd - panel->cny;
+
+ *fs = panel->xfs*plx + panel->yfs*ply;
+ *ss = panel->xss*plx + panel->yss*ply;
+
+ *fs += panel->min_fs;
+ *ss += panel->min_ss;
+
+ /* Now, is this on this panel? */
+ if ( *fs < panel->min_fs ) return 3;
+ if ( *fs > panel->max_fs ) return 3;
+ if ( *ss < panel->min_ss ) return 3;
+ if ( *ss > panel->max_ss ) return 3;
+
+ return 0;
+
+}
+
+//FIXME: should this function be in detector.h?
+struct rvec nearest_bragg(struct image image, struct rvec q)
+{
+
+ struct rvec g;
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ double hd, kd, ld;
+ double h, k, l;
+
+ /* miller indices of nearest Bragg reflection */
+ cell_get_cartesian(image.indexed_cell, &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz);
+
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+
+ h = lrint(hd);
+ k = lrint(kd);
+ l = lrint(ld);
+
+ /* now get scattering vector for reflectin [hkl]
+ * this means solving the equation U*q = h */
+ double U[] = {ax, ay, az,
+ bx, by, bz,
+ cx, cy, cz};
+
+ double hvec[] = {h,k,l};
+
+ gsl_matrix_view m
+ = gsl_matrix_view_array (U, 3, 3);
+
+ gsl_vector_view b
+ = gsl_vector_view_array (hvec, 3);
+
+ gsl_vector *x = gsl_vector_alloc (3);
+
+ int s;
+
+ gsl_permutation * perm = gsl_permutation_alloc (3);
+ gsl_linalg_LU_decomp (&m.matrix, perm, &s);
+ gsl_linalg_LU_solve (&m.matrix, perm, &b.vector, x);
+
+ /* outgoing wavevector for [hkl] */
+ g.u = x->data[0];
+ g.v = x->data[1];
+ g.w = x->data[2];
+
+ return g;
+
+}
int main(int argc, char *argv[])
@@ -214,13 +311,9 @@ int main(int argc, char *argv[])
struct panel * p;
struct imagefeature * thisFeature;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- double hd, kd, ld; /* Indices with decimal places */
- signed int h, k, l;
struct rvec q;
double twotheta;
+ struct rvec g;
double fs, ss; /* observed peaks */
double pfs, pss; /* predicted peaks */
double dfs, dss; /* observed - predicted */
@@ -243,90 +336,19 @@ int main(int argc, char *argv[])
}
if ( p->no_index ) continue;
-
-
/* now determine the predicted peak position */
/* scattering vector of this peak */
q = get_q(&image, fs, ss, &twotheta, 1.0/image.lambda);
- //fail = calculate_projected_peak();
+ /* scattering vector of nearest bragg peak */
+ g = nearest_bragg(image, q);
- /* miller indices of nearest Bragg reflection */
- cell_get_cartesian(image.indexed_cell, &ax, &ay, &az,
- &bx, &by, &bz,
- &cx, &cy, &cz);
-
- hd = q.u * ax + q.v * ay + q.w * az;
- kd = q.u * bx + q.v * by + q.w * bz;
- ld = q.u * cx + q.v * cy + q.w * cz;
-
- h = lrint(hd);
- k = lrint(kd);
- l = lrint(ld);
-
- /* now get scattering vector for reflectin [hkl]
- * this means solving the equation U*q = h */
-
- double U[] = {ax, ay, az,
- bx, by, bz,
- cx, cy, cz};
-
- double hvec[] = {h,k,l};
-
- gsl_matrix_view m
- = gsl_matrix_view_array (U, 3, 3);
-
- gsl_vector_view b
- = gsl_vector_view_array (hvec, 3);
-
- gsl_vector *x = gsl_vector_alloc (3);
-
- int s;
-
- gsl_permutation * perm = gsl_permutation_alloc (3);
- gsl_linalg_LU_decomp (&m.matrix, perm, &s);
- gsl_linalg_LU_solve (&m.matrix, perm, &b.vector, x);
-
-
- /* outgoing wavevector for [hkl] */
- double gx = x->data[0];
- double gy = x->data[1];
- double gz = x->data[2];
-
- double kk;
- double xd, yd, cl;
- double plx, ply;
-
- kk = 1/image.lambda;
- const double den = kk + gz;
-
- /* Camera length for this panel */
- cl = p->clen;
-
- /* Coordinates of peak relative to central beam, in m */
- xd = cl * gx / den;
- yd = cl * gy / den;
-
- /* Convert to pixels */
- xd *= p->res;
- yd *= p->res;
-
- /* Convert to relative to the panel corner */
- plx = xd - p->cnx;
- ply = yd - p->cny;
-
- pfs = p->xfs*plx + p->yfs*ply;
- pss = p->xss*plx + p->yss*ply;
-
- pfs += p->min_fs;
- pss += p->min_ss;
-
- /* Now, is this on this panel? */
- if ( fs < p->min_fs ) continue;
- if ( fs > p->max_fs ) continue;
- if ( ss < p->min_ss ) continue;
- if ( ss > p->max_ss ) continue;
+ /* coordinate of this predicted peak */
+ fail = calculate_projected_peak(p,g,1/image.lambda,&pfs,&pss);
+
+ /* check for error, e.g. out of panel */
+ if ( fail != 0 ) continue;
/* Finally, we have the shift in position of this peak */
dfs = pfs - fs;
@@ -344,16 +366,13 @@ int main(int argc, char *argv[])
} /* end loop through stream chunks */
-
-
/* calculate weighted average shift in peak positions */
for (pi=0; pi < image.det->n_panels; pi++) {
meanShiftFS[pi] = weightedSumFS[pi]/summedWeights[pi];
meanShiftSS[pi] = weightedSumSS[pi]/summedWeights[pi];
}
- /* now generate a new geometry file */
- /* first populate the image structure with new geometry */
+ /* first populate the image structure with new geometry info */
for (pi=0; pi < image.det->n_panels; pi++) {
p = image.det->panels[pi];
@@ -363,8 +382,6 @@ int main(int argc, char *argv[])
if ( peaksFound[pi] >= minpeaks ) {
- //printf("meanShift: %f %f\n",meanShiftFS[pi],meanShiftSS[pi]);
-
/* convert shifts from raw coords to lab frame */
xsh = meanShiftFS[pi]*p.fsx + meanShiftSS[pi]*p.ssx;
ysh = meanShiftFS[pi]*p.fsy + meanShiftSS[pi]*p.ssy;
@@ -373,30 +390,21 @@ int main(int argc, char *argv[])
cnx = p.cnx + xsh;
cny = p.cny + ysh;
- //printf("new panel shift: %f %f\n",cnx,cny);
-
-
} else {
-
- /* not refined: use original coordinates */
+ /* not refined? use original coordinates */
cnx = p.cnx;
cny = p.cny;
-
}
- /* now dump this refined panel goemetry into a file */
-
image.det->panels[pi].cnx = cnx;
image.det->panels[pi].cny = cny;
if ( peaksFound[pi] < minpeaks) image.det->panels[pi].no_index = 1;
printf("panel %s, # peaks: %10d, mean shifts: %f %f\n",p.name, peaksFound[pi],xsh,ysh);
- //fprintf_panel(outfh,&p);
- //fprintf(outfh,"\n");
-
}
+ /* write the new geometry file */
print_detector_geometry(outfh,&image);
} else {
@@ -406,11 +414,8 @@ int main(int argc, char *argv[])
}
- if ( !(outputfile == NULL) ) {
- fclose(outfh);
- }
-
+ fclose(outfh);
printf("Done!\n");
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 02:24:22 +0000