From c080e2c362628da1a68ad17c117f303945724006 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 28 Feb 2011 17:11:00 +0100 Subject: First round of fallout from new geometry format --- src/detector.c | 60 ++++++++++++++++++++++++++++++++++++----------------- src/detector.h | 9 ++++++-- src/displaywindow.c | 4 +++- src/geometry.c | 41 ++++++++++++++++++++++-------------- src/mosflm.c | 18 ++++++++++------ src/peaks.c | 8 +++++-- 6 files changed, 94 insertions(+), 46 deletions(-) (limited to 'src') 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; in_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; iimage->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; pn_panels; p++ ) { + for ( i=0; in_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; ifeatures, 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); -- cgit v1.2.3