Add the possibility to use a circular (instead of 'crystallographic') integration domain

5 files changed, 37 insertions, 15 deletions

diff --git a/src/indexamajig.c b/src/indexamajig.c
index 39905fb4..281294f3 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -362,7 +362,7 @@ static struct process_result process_image(struct process_args *pargs)
 
 	/* Sanity check */
 	if ( pargs->config_sanity
-	  && !peak_sanity_check(&image, image.indexed_cell) ) {
+	  && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) {
 		STATUS("Failed peak sanity check.\n");
 		goto done;
 	} else {
@@ -373,7 +373,8 @@ static struct process_result process_image(struct process_args *pargs)
 	if ( config_nearbragg ) {
 		output_intensities(&image, image.indexed_cell,
 		                   pargs->output_mutex, config_polar,
-		                   pargs->config_sa, pargs->config_closer);
+		                   pargs->config_sa, pargs->config_closer,
+		                   0, 0.1);
 	}
 
 	simage = get_simage(&image, config_alternate);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 4d4ba10f..0d414653 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -417,8 +417,8 @@ int main(int argc, char *argv[])
 		record_image(&image, !config_nonoise);
 
 		if ( config_nearbragg ) {
-			find_projected_peaks(&image, cell);
-			output_intensities(&image, cell, NULL, 0, 1, 0);
+			find_projected_peaks(&image, cell, 0, 0.1);
+			output_intensities(&image, cell, NULL, 0, 1, 0, 0, 0.1);
 		}
 
 		if ( powder_fn != NULL ) {
diff --git a/src/peaks.c b/src/peaks.c
index fdf5b961..5cfdbdde 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -451,7 +451,8 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
 }
 
 
-int find_projected_peaks(struct image *image, UnitCell *cell)
+int find_projected_peaks(struct image *image, UnitCell *cell,
+                         int circular_domain, double domain_r)
 {
 	int x, y;
 	double ax, ay, az;
@@ -459,11 +460,15 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 	double cx, cy, cz;
 	struct reflhit *hits;
 	int n_hits = 0;
+	double alen, blen, clen;
 
 	hits = malloc(sizeof(struct reflhit)*MAX_HITS);
 	if ( hits == NULL ) return 0;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	alen = modulus(ax, ay, az);
+	blen = modulus(bx, by, bz);
+	clen = modulus(cx, cy, cz);
 
 	fesetround(1);  /* Round towards nearest */
 	for ( x=0; x<image->width; x++ ) {
@@ -490,8 +495,17 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 		dh = hd - h;
 		dk = kd - k;
 		dl = ld - l;
-		dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
-		if ( dist > 0.1 ) continue;
+
+		if ( circular_domain ) {
+			/* Circular integration domain */
+			dist = sqrt(pow(dh*alen, 2.0) + pow(dk*blen, 2.0)
+			                              + pow(dl*clen, 2.0));
+			if ( dist > domain_r ) continue;
+		} else {
+			/* "Crystallographic" integration domain */
+			dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
+			if ( dist > domain_r ) continue;
+		}
 
 		for ( j=0; j<n_hits; j++ ) {
 			if ( (hits[j].h == h) && (hits[j].k == k)
@@ -530,12 +544,13 @@ int find_projected_peaks(struct image *image, UnitCell *cell)
 }
 
 
-int peak_sanity_check(struct image *image, UnitCell *cell)
+int peak_sanity_check(struct image *image, UnitCell *cell,
+                      int circular_domain, double domain_r)
 {
 	int i;
 	int n_sane = 0;
 
-	find_projected_peaks(image, cell);
+	find_projected_peaks(image, cell, circular_domain, domain_r);
 	if ( image->n_hits == 0 ) return 0;  /* Failed sanity check: no peaks */
 
 	for ( i=0; i<image->n_hits; i++ ) {
@@ -564,7 +579,7 @@ int peak_sanity_check(struct image *image, UnitCell *cell)
 
 void output_intensities(struct image *image, UnitCell *cell,
                         pthread_mutex_t *mutex, int polar, int sa,
-                        int use_closer)
+                        int use_closer, int circular_domain, double domain_r)
 {
 	int i;
 	int n_found;
@@ -577,7 +592,9 @@ void output_intensities(struct image *image, UnitCell *cell,
 	double csx, csy, csz;
 	double a, b, c, al, be, ga;
 
-	if ( image->n_hits == 0 ) find_projected_peaks(image, cell);
+	if ( image->n_hits == 0 ) {
+		find_projected_peaks(image, cell, circular_domain, domain_r);
+	}
 	if ( image->n_hits == 0 ) return;
 
 	/* Get exclusive access to the output stream if necessary */
diff --git a/src/peaks.h b/src/peaks.h
index 68e104a9..2d8c1db5 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -23,9 +23,12 @@ extern void search_peaks(struct image *image, float threshold);
 extern void dump_peaks(struct image *image, pthread_mutex_t *mutex);
 extern void output_intensities(struct image *image, UnitCell *cell,
                                pthread_mutex_t *mutex, int polar, int sa,
-                               int use_closer);
-extern int peak_sanity_check(struct image *image, UnitCell *cell);
-extern int find_projected_peaks(struct image *image, UnitCell *cell);
+                               int use_closer, int circular_domain,
+                               double domain_r);
+extern int peak_sanity_check(struct image *image, UnitCell *cell,
+                             int circular_domain, double domain_r);
+extern int find_projected_peaks(struct image *image, UnitCell *cell,
+                                int circular_domain, double domain_r);
 extern int integrate_peak(struct image *image, int xp, int yp,
                           float *xc, float *yc, float *intensity,
                           int do_polar, int do_sa);
diff --git a/src/reintegrate.c b/src/reintegrate.c
index 3cb94730..3a94c48f 100644
--- a/src/reintegrate.c
+++ b/src/reintegrate.c
@@ -103,7 +103,8 @@ static void process_image(struct process_args *pargs)
 
 	output_intensities(&image, pargs->cell,
 		                   pargs->output_mutex, 1,
-		                   1, 0);
+		                   1, 0,
+		                   0, 0.1);
 
 	free(image.data);
 	cell_free(pargs->cell);