First round of fallout from new geometry format

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);