aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-02-28 17:11:00 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:16 +0100
commitc080e2c362628da1a68ad17c117f303945724006 (patch)
treeb24e0c013e34dc2d339d91371ce0c3655a037831
parent0d2fd58266f87fb20aceb017c68b0e455ab4baee (diff)
First round of fallout from new geometry format
-rw-r--r--src/detector.c60
-rw-r--r--src/detector.h9
-rw-r--r--src/displaywindow.c4
-rw-r--r--src/geometry.c41
-rw-r--r--src/mosflm.c18
-rw-r--r--src/peaks.c8
6 files changed, 94 insertions, 46 deletions
diff --git a/src/detector.c b/src/detector.c
index 7d2cb513..703a7189 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -82,8 +82,8 @@ struct rvec get_q(struct image *image, double fs, double ss,
xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx;
ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy;
- rx = (xs + p->cx) / p->res;
- ry = (ys + p->cy) / p->res;
+ rx = (xs + p->cnx) / p->res;
+ ry = (ys + p->cny) / p->res;
/* Calculate q-vector for this sub-pixel */
r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
@@ -113,11 +113,8 @@ double get_tt(struct image *image, double fs, double ss)
xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx;
ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy;
- rx = (xs + p->cx) / p->res;
- ry = (ys + p->cy) / p->res;
-
- /* Calculate q-vector for this sub-pixel */
- r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+ rx = (xs + p->cnx) / p->res;
+ ry = (ys + p->cny) / p->res;
r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
@@ -150,6 +147,7 @@ void record_image(struct image *image, int do_poisson)
double cf;
double intensity, sa;
double pix_area, Lsq;
+ double xs, ys, rx, ry;
double dsq, proj_area;
struct panel *p;
@@ -174,8 +172,11 @@ void record_image(struct image *image, int do_poisson)
proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
/* Calculate distance from crystal to pixel */
- dsq = pow(((double)x - p->cx) / p->res, 2.0);
- dsq += pow(((double)y - p->cy) / p->res, 2.0);
+ xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
+ ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
+ rx = (xs + p->cnx) / p->res;
+ ry = (ys + p->cny) / p->res;
+ dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
/* Projected area of pixel divided by distance squared */
sa = proj_area / (dsq + Lsq);
@@ -324,8 +325,8 @@ struct detector *get_detector_geometry(const char *filename)
det->panels[i].min_ss = -1;
det->panels[i].max_fs = -1;
det->panels[i].max_ss = -1;
- det->panels[i].cx = -1;
- det->panels[i].cy = -1;
+ det->panels[i].cnx = -1;
+ det->panels[i].cny = -1;
det->panels[i].clen = -1;
det->panels[i].res = -1;
det->panels[i].badrow = '-';
@@ -374,9 +375,9 @@ struct detector *get_detector_geometry(const char *filename)
} else if ( strcmp(path[1], "max_ss") == 0 ) {
det->panels[np].max_ss = atof(bits[2]);
} else if ( strcmp(path[1], "corner_x") == 0 ) {
- det->panels[np].cx = atof(bits[2]);
+ det->panels[np].cnx = atof(bits[2]);
} else if ( strcmp(path[1], "corner_y") == 0 ) {
- det->panels[np].cy = atof(bits[2]);
+ det->panels[np].cny = atof(bits[2]);
} else if ( strcmp(path[1], "clen") == 0 ) {
char *end;
@@ -462,12 +463,12 @@ struct detector *get_detector_geometry(const char *filename)
" panel %i\n", i);
reject = 1;
}
- if ( det->panels[i].cx == -1 ) {
+ if ( det->panels[i].cnx == -1 ) {
ERROR("Please specify the corner X coordinate for"
" panel %i\n", i);
reject = 1;
}
- if ( det->panels[i].cy == -1 ) {
+ if ( det->panels[i].cny == -1 ) {
ERROR("Please specify the corner Y coordinate for"
" panel %i\n", i);
reject = 1;
@@ -509,6 +510,27 @@ out:
det->max_fs = max_fs;
det->max_ss = max_ss;
+ /* Calculate matrix inverse */
+ for ( i=0; i<det->n_panels; i++ ) {
+
+ struct panel *p;
+ double d;
+
+ p = &det->panels[i];
+
+ if ( p->fsx*p->ssy == p->ssx*p->fsy ) {
+ ERROR("Panel %i transformation singular.\n", i);
+ reject = 1;
+ }
+
+ d = (double)p->fsx*p->ssy - p->ssx*p->fsy;
+ p->xfs = p->ssy / d;
+ p->yfs = -p->ssx / d;
+ p->xss = -p->fsy / d;
+ p->yss = p->fsx / d;
+
+ }
+
if ( reject ) return NULL;
fclose(fh);
@@ -537,8 +559,8 @@ struct detector *simple_geometry(const struct image *image)
geom->panels[0].max_fs = image->width-1;
geom->panels[0].min_ss = 0;
geom->panels[0].max_ss = image->height-1;
- geom->panels[0].cx = -image->width / 2.0;
- geom->panels[0].cy = -image->height / 2.0;
+ geom->panels[0].cnx = -image->width / 2.0;
+ geom->panels[0].cny = -image->height / 2.0;
geom->panels[0].fsx = 1;
geom->panels[0].fsy = 0;
geom->panels[0].ssx = 0;
@@ -556,8 +578,8 @@ static void check_extents(struct panel p, double *min_x, double *min_y,
xs = fs*p.fsx + ss*p.ssx;
ys = fs*p.fsy + ss*p.ssy;
- rx = xs + p.cx;
- ry = ys + p.cy;
+ rx = xs + p.cnx;
+ ry = ys + p.cny;
if ( rx > *max_x ) *max_x = rx;
if ( ry > *max_y ) *max_y = ry;
diff --git a/src/detector.h b/src/detector.h
index e9b2112f..a232b712 100644
--- a/src/detector.h
+++ b/src/detector.h
@@ -29,8 +29,8 @@ struct panel
int max_fs; /* Largest FS value considered to be in this panel */
int min_ss; /* ... and so on */
int max_ss;
- double cx; /* Location of corner (min_fs,min_ss) in pixels */
- double cy;
+ double cnx; /* Location of corner (min_fs,min_ss) in pixels */
+ double cny;
double clen; /* Camera length in metres */
char *clen_from;
double res; /* Resolution in pixels per metre */
@@ -42,6 +42,11 @@ struct panel
signed int fsy;
signed int ssx;
signed int ssy;
+
+ double xfs;
+ double yfs;
+ double xss;
+ double yss;
};
struct detector
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 51a5cb9d..b542f269 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -202,7 +202,8 @@ static gboolean displaywindow_expose(GtkWidget *da, GdkEventExpose *event,
cairo_set_matrix(cr, &basic_m);
/* Move to the right location */
- cairo_translate(cr, p.cx/dw->binning, p.cy/dw->binning);
+ cairo_translate(cr, p.cnx/dw->binning,
+ p.cny/dw->binning);
/* Twiddle directions according to matrix */
cairo_matrix_init(&m, p.fsx, p.fsy, p.ssx, p.ssy,
@@ -248,6 +249,7 @@ static gboolean displaywindow_expose(GtkWidget *da, GdkEventExpose *event,
return FALSE;
}
+ cairo_identity_matrix(cr);
for ( i=0; i<image_feature_count(dw->image->features); i++ ) {
double x, y;
diff --git a/src/geometry.c b/src/geometry.c
index 33527f71..f8b74148 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -31,45 +31,54 @@
static signed int locate_peak(double x, double y, double z, double k,
struct detector *det, double *xdap, double *ydap)
{
- int p;
+ int i;
signed int found = -1;
const double den = k + z;
*xdap = -1; *ydap = -1;
- for ( p=0; p<det->n_panels; p++ ) {
+ for ( i=0; i<det->n_panels; i++ ) {
double xd, yd, cl;
- double xda, yda;
+ double fs, ss, plx, ply;
+ struct panel *p;
+
+ p = &det->panels[i];
/* Camera length for this panel */
- cl = det->panels[p].clen;
+ cl = p->clen;
/* Coordinates of peak relative to central beam, in m */
xd = cl * x / den;
yd = cl * y / den;
/* Convert to pixels */
- xd *= det->panels[p].res;
- yd *= det->panels[p].res;
+ xd *= p->res;
+ yd *= p->res;
+
+ /* Convert to relative to the panel corner */
+ plx = xd - p->cnx;
+ ply = yd - p->cny;
+
+ fs = p->xfs*plx + p->yfs*ply;
+ ss = p->xss*plx + p->yss*ply;
- /* Add the coordinates of the central beam */
- xda = xd + det->panels[p].cx;
- yda = yd + det->panels[p].cy;
+ fs += p->min_fs;
+ ss += p->min_ss;
/* Now, is this on this panel? */
- if ( xda < det->panels[p].min_fs ) continue;
- if ( xda > det->panels[p].max_fs ) continue;
- if ( yda < det->panels[p].min_ss ) continue;
- if ( yda > det->panels[p].max_ss ) continue;
+ if ( fs < p->min_fs ) continue;
+ if ( fs > p->max_fs ) continue;
+ if ( ss < p->min_ss ) continue;
+ if ( ss > p->max_ss ) continue;
/* If peak appears on multiple panels, reject it */
if ( found != -1 ) return -1;
/* Woohoo! */
- found = p;
- *xdap = xda;
- *ydap = yda;
+ found = i;
+ *xdap = fs;
+ *ydap = ss;
}
diff --git a/src/mosflm.c b/src/mosflm.c
index bfaca2fb..562b58d7 100644
--- a/src/mosflm.c
+++ b/src/mosflm.c
@@ -199,19 +199,25 @@ static void write_spt(struct image *image, const char *filename)
for ( i=0; i<nPeaks; i++ ) {
struct imagefeature *f;
+ struct panel *p;
+ double xs, ys, rx, ry;
f = image_get_feature(image->features, i);
if ( f == NULL ) continue;
- struct panel *pan;
- pan = find_panel(image->det, f->x, f->y);
- if ( pan == NULL ) continue;
+ p = find_panel(image->det, f->x, f->y);
+ if ( p == NULL ) continue;
- pix = 1000/pan->res; /* pixel size in mm */
+ pix = 1000.0/p->res; /* pixel size in mm */
height = f->intensity;
- sptlines[i].x = (f->y - pan->cy)*pix*fclen/pan->clen/1000;
- sptlines[i].y = -(f->x - pan->cx)*pix*fclen/pan->clen/1000;
+ xs = (f->x-p->min_fs)*p->fsx + (f->y-p->min_ss)*p->ssx;
+ ys = (f->x-p->min_fs)*p->fsy + (f->y-p->min_ss)*p->ssy;
+ rx = xs + p->cnx;
+ ry = ys + p->cny;
+
+ sptlines[i].x = ry*pix*fclen/p->clen/1000.0;
+ sptlines[i].y = -rx*pix*fclen/p->clen/1000.0;
sptlines[i].h = height;
sptlines[i].s = sigma;
diff --git a/src/peaks.c b/src/peaks.c
index 676a4387..640188e7 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -857,6 +857,7 @@ void output_pixels(struct image *image, UnitCell *cell,
double pix_area, Lsq, proj_area, dsq, sa;
double phi, pa, pb, pol;
float tt = 0.0;
+ double xs, ys, rx, ry;
/* Veto if we want to integrate a bad region */
if ( image->flags != NULL ) {
@@ -881,8 +882,11 @@ void output_pixels(struct image *image, UnitCell *cell,
proj_area = pix_area * cos(tt);
/* Calculate distance from crystal to pixel */
- dsq = pow(((double)x - p->cx) / p->res, 2.0);
- dsq += pow(((double)y - p->cy) / p->res, 2.0);
+ xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
+ ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
+ rx = (xs + p->cnx) / p->res;
+ ry = (ys + p->cny) / p->res;
+ dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
/* Projected area of pixel / distance squared */
sa = 1.0e7 * proj_area / (dsq + Lsq);