render_hkl: Add alternative weighting schemes

2 files changed, 94 insertions, 49 deletions

diff --git a/src/reflections.c b/src/reflections.c
index 8a9ad163..2711936b 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -115,15 +115,15 @@ double *read_reflections(const char *filename, unsigned int *counts,
 
 		char line[1024];
 		signed int h, k, l;
-		float intensity, ph, res;
+		float intensity, ph, res, sigma;
 		char phs[1024];
 		int r;
 		int cts;
 
 		rval = fgets(line, 1023, fh);
-		r = sscanf(line, "%i %i %i %f %s %f %i",
-		           &h, &k, &l, &intensity, phs, &res, &cts);
-		if ( r != 5 ) continue;
+		r = sscanf(line, "%i %i %i %f %s %f %f %i",
+		           &h, &k, &l, &intensity, phs, &sigma, &res, &cts);
+		if ( r != 8 ) continue;
 
 		set_intensity(ref, h, k, l, intensity);
 		if ( phases != NULL ) {
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 7bfe97b3..ff955084 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -28,6 +28,11 @@
 #include "povray.h"
 #include "symmetry.h"
 
+enum {
+	WGHT_I,
+	WGHT_SQRTI,
+	WGHT_COUNTS,
+};
 
 static void show_help(const char *s)
 {
@@ -35,25 +40,28 @@ static void show_help(const char *s)
 	printf(
 "Render intensity lists in various ways.\n"
 "\n"
-"  -h, --help            Display this help message.\n"
-"      --povray          Render a 3D animation using POV-ray.\n"
-"      --zone-axis       Render a 2D zone axis pattern.\n"
-"      --boost=<val>     Squash colour scale by <val>.\n"
-"  -j <n>                Run <n> instances of POV-ray in parallel.\n"
-"  -p, --pdb=<file>      PDB file from which to get the unit cell.\n"
-"  -y, --symmetry=<sym>  Expand reflections according to point group <sym>.\n"
-"      --sqrt            Plot square roots of intensities, rather than\n"
-"                         the intensities themselves.\n"
+"  -h, --help              Display this help message.\n"
+"      --povray            Render a 3D animation using POV-ray.\n"
+"      --zone-axis         Render a 2D zone axis pattern.\n"
+"      --boost=<val>       Squash colour scale by <val>.\n"
+"  -j <n>                  Run <n> instances of POV-ray in parallel.\n"
+"  -p, --pdb=<file>        PDB file from which to get the unit cell.\n"
+"  -y, --symmetry=<sym>    Expand reflections according to point group <sym>.\n"
+"  -w  --weighting=<wght>  Colour/shade the reciprocal lattice points\n"
+"                           according to:\n"
+"                            I     : the intensity of the reflection.\n"
+"                            sqrtI : the square root of the intensity.\n"
+"                            count : the number of counts for the reflection.\n"
 );
 }
 
 
 static void render_za(UnitCell *cell, double *ref, unsigned int *c,
-                      double boost, const char *sym, int config_sqrt)
+                      double boost, const char *sym, int wght)
 {
 	cairo_surface_t *surface;
 	cairo_t *dctx;
-	double max_u, max_v, max_res, max_intensity;
+	double max_u, max_v, max_res, max_val;
 	double scale_u, scale_v, scale;
 	double sep_u, sep_v, max_r;
 	double u, v;
@@ -83,7 +91,7 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
 	cairo_set_source_rgb(dctx, 0.0, 0.0, 0.0);
 	cairo_fill(dctx);
 
-	max_u = 0.0;  max_v = 0.0;  max_intensity = 0.0;
+	max_u = 0.0;  max_v = 0.0;  max_val = 0.0;
 	max_res = 0.0;  max_h = 0;  max_k = 0;
 
 	/* Work out reciprocal lattice spacings and angles for this cut */
@@ -101,38 +109,44 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
 	for ( h=-INDMAX; h<INDMAX; h++ ) {
 	for ( k=-INDMAX; k<INDMAX; k++ ) {
 
-		double u, v, intensity, res;
+		double u, v, val, res;
 		int ct;
 		int nequiv, p;
 
 		ct = lookup_count(c, h, k, 0);
 		if ( ct < 1 ) continue;
 
-		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
-		if ( config_sqrt ) {
-			intensity = (intensity>0.0) ? sqrt(intensity) : 0.0;
+		switch ( wght ) {
+		case WGHT_I :
+			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			break;
+		case WGHT_SQRTI :
+			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			val = (val>0.0) ? sqrt(val) : 0.0;
+			break;
+		case WGHT_COUNTS :
+			val = (float)ct;
+			break;
 		}
-		if ( intensity != 0 ) {
-
-			nequiv = num_equivs(h, k, 0, sym);
-			for ( p=0; p<nequiv; p++ ) {
-
-				signed int he, ke, le;
-				get_equiv(h, k, 0, &he, &ke, &le, sym, p);
-
-				u =  (double)he*as*sin(theta);
-				v =  (double)he*as*cos(theta) + ke*bs;
-				if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
-				if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
-				if ( fabs(intensity) > fabs(max_intensity) ) {
-					max_intensity = fabs(intensity);
-				}
-				if ( fabs(he) > max_h ) max_h = fabs(he);
-				if ( fabs(ke) > max_k ) max_k = fabs(ke);
-				res = resolution(cell, he, ke, 0);
-				if ( res > max_res ) max_res = res;
 
+		nequiv = num_equivs(h, k, 0, sym);
+		for ( p=0; p<nequiv; p++ ) {
+
+			signed int he, ke, le;
+			get_equiv(h, k, 0, &he, &ke, &le, sym, p);
+
+			u =  (double)he*as*sin(theta);
+			v =  (double)he*as*cos(theta) + ke*bs;
+			if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
+			if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
+			if ( fabs(val) > fabs(max_val) ) {
+				max_val = fabs(val);
 			}
+			if ( fabs(he) > max_h ) max_h = fabs(he);
+			if ( fabs(ke) > max_k ) max_k = fabs(ke);
+			res = resolution(cell, he, ke, 0);
+			if ( res > max_res ) max_res = res;
+
 		}
 
 	}
@@ -142,7 +156,7 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
 	printf("Maximum resolution is 1/d = %5.3f nm^-1, d = %5.3f nm\n",
 	       max_res*2.0, 1.0/(max_res*2.0));
 
-	if ( max_intensity <= 0.0 ) {
+	if ( max_val <= 0.0 ) {
 		max_r = 4.0;
 		goto out;
 	}
@@ -160,18 +174,27 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
 	for ( h=-INDMAX; h<INDMAX; h++ ) {
 	for ( k=-INDMAX; k<INDMAX; k++ ) {
 
-		double u, v, intensity, val;
+		double u, v, val;
 		int ct;
 		int nequiv, p;
 
 		ct = lookup_count(c, h, k, 0);
-		if ( ct < 1 ) continue;
+		if ( ct < 1 ) continue;  /* Must have at least one count */
 
-		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
-		if ( config_sqrt ) {
-			intensity = (intensity>0.0) ? sqrt(intensity) : 0.0;
+		switch ( wght ) {
+		case WGHT_I :
+			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			break;
+		case WGHT_SQRTI :
+			val = lookup_intensity(ref, h, k, 0) / (float)ct;
+			val = (val>0.0) ? sqrt(val) : 0.0;
+			break;
+		case WGHT_COUNTS :
+			val = (float)ct;
+			break;
 		}
-		val = boost*intensity/max_intensity;
+
+		val = boost*val/max_val;
 
 		nequiv = num_equivs(h, k, 0, sym);
 		for ( p=0; p<nequiv; p++ ) {
@@ -251,6 +274,8 @@ int main(int argc, char *argv[])
 	int r = 0;
 	double boost = 1.0;
 	char *sym = NULL;
+	char *weighting = NULL;
+	int wght;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -260,12 +285,13 @@ int main(int argc, char *argv[])
 		{"pdb",                1, NULL,               'p'},
 		{"boost",              1, NULL,               'b'},
 		{"symmetry",           1, NULL,               'y'},
-		{"sqrt",               0, &config_sqrt,        1},
+		{"weighting",          1, NULL,               'w'},
+		{"counts",             0, &config_sqrt,        1},
 		{0, 0, NULL, 0}
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hj:p:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "hj:p:w:", longopts, NULL)) != -1) {
 
 		switch (c) {
 		case 'h' :
@@ -288,6 +314,10 @@ int main(int argc, char *argv[])
 			sym = strdup(optarg);
 			break;
 
+		case 'w' :
+			weighting = strdup(optarg);
+			break;
+
 		case 0 :
 			break;
 
@@ -305,6 +335,21 @@ int main(int argc, char *argv[])
 		sym = strdup("1");
 	}
 
+	if ( weighting == NULL ) {
+		weighting = strdup("I");
+	}
+
+	if ( strcmp(weighting, "I") == 0 ) {
+		wght = WGHT_I;
+	} else if ( strcmp(weighting, "sqrtI") == 0 ) {
+		wght = WGHT_SQRTI;
+	} else if ( strcmp(weighting, "count") == 0 ) {
+		wght = WGHT_COUNTS;
+	} else {
+		ERROR("Unrecognised weighting '%s'\n", weighting);
+		return 1;
+	}
+
 	infile = argv[optind];
 
 	cell = load_cell_from_pdb(pdb);
@@ -322,7 +367,7 @@ int main(int argc, char *argv[])
 	if ( config_povray ) {
 		r = povray_render_animation(cell, ref, cts, nproc);
 	} else if ( config_zoneaxis ) {
-		render_za(cell, ref, cts, boost, sym, config_sqrt);
+		render_za(cell, ref, cts, boost, sym, wght);
 	} else {
 		ERROR("Try again with either --povray or --zone-axis.\n");
 	}