partialator: Use all reflections (not just the ones found the first time around)

doc/reference/CrystFEL-sections.txt

@@ -6,8 +6,11 @@ RefListIterator
 <SUBSECTION>
 reflist_new
 reflist_free
+reflection_new
+reflection_free
 <SUBSECTION>
 add_refl
+add_refl_to_list
 <SUBSECTION>
 first_refl
 next_refl

src/geometry.c

@@ -122,7 +122,7 @@ static double partiality(double r1, double r2, double r)
 }
 
 
-static Reflection *check_reflection(struct image *image, RefList *reflections,
+static Reflection *check_reflection(struct image *image,
                                     signed int h, signed int k, signed int l,
                                     double asx, double asy, double asz,
                                     double bsx, double bsy, double bsz,
@@ -216,7 +216,7 @@ static Reflection *check_reflection(struct image *image, RefList *reflections,
 	if ( p == -1 ) return NULL;
 
 	/* Add peak to list */
-	refl = add_refl(reflections, h, k, l);
+	refl = reflection_new(h, k, l);
 	set_detector_pos(refl, 0.0, xda, yda);
 	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
 	set_symmetric_indices(refl, h, k, l);
@@ -274,8 +274,16 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 	for ( h=-hmax; h<=hmax; h++ ) {
 	for ( k=-kmax; k<=kmax; k++ ) {
 	for ( l=-lmax; l<=lmax; l++ ) {
-		check_reflection(image, reflections, h, k, l,
-		                 asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+		Reflection *refl;
+
+		refl = check_reflection(image, h, k, l,
+		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
+
+		if ( refl != NULL ) {
+			refl = add_refl_to_list(refl, reflections);
+		}
+
 	}
 	}
 	}
@@ -285,8 +293,7 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 
 
 /* Calculate partialities and apply them to the image's reflections */
-void update_partialities(struct image *image,
-                         int *n_expected, int *n_found, int *n_notfound)
+void update_partialities(struct image *image)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -295,10 +302,6 @@ void update_partialities(struct image *image,
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 
-	if ( n_expected != NULL ) *n_expected = 0;
-	if ( n_found != NULL) *n_found = 0;
-	if ( n_notfound != NULL) *n_notfound = 0;
-
 	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
 	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
@@ -314,19 +317,15 @@ void update_partialities(struct image *image,
 		signed int h, k, l;
 		int clamp1, clamp2;
 
-		if ( n_expected != NULL ) (*n_expected)++;
-
 		get_symmetric_indices(refl, &h, &k, &l);
 
-		vals = check_reflection(image, predicted, h, k, l,
+		vals = check_reflection(image, h, k, l,
 		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
 
 		if ( vals == NULL ) {
-			if ( n_notfound != NULL) (*n_notfound)++;
 			set_redundancy(refl, 0);
 			continue;
 		}
-		if ( n_found != NULL) (*n_found)++;
 		set_redundancy(refl, 1);
 
 		/* Transfer partiality stuff */

src/geometry.h

@@ -21,7 +21,6 @@
 
 extern RefList *find_intersections(struct image *image, UnitCell *cell);
 
-extern void update_partialities(struct image *image,
-                                int *n_expected, int *n_found, int *n_notfound);
+extern void update_partialities(struct image *image);
 
 #endif	/* GEOMETRY_H */

src/partialator.c

@@ -269,9 +269,6 @@ int main(int argc, char *argv[])
 	int n_iter = 10;
 	struct beam_params *beam = NULL;
 	RefList *full;
-	int n_found = 0;
-	int n_expected = 0;
-	int n_notfound = 0;
 	int n_usable_patterns = 0;
 	int nobs;
 	char *reference_file = NULL;
@@ -435,7 +432,6 @@ int main(int argc, char *argv[])
 
 		RefList *as;
 		struct image *cur;
-		int nn_expected, nn_found, nn_notfound;
 
 		cur = &images[n_usable_patterns];
 
@@ -474,10 +470,7 @@ int main(int argc, char *argv[])
 		reflist_free(cur->reflections);
 		cur->reflections = as;
 
-		update_partialities(cur, &nn_expected, &nn_found, &nn_notfound);
-		n_expected += nn_expected;
-		n_found += nn_found;
-		n_notfound += nn_notfound;
+		update_partialities(cur);
 
 		nobs += select_scalable_reflections(cur->reflections,
 		                                    reference);
@@ -487,8 +480,6 @@ int main(int argc, char *argv[])
 
 	}
 	fclose(fh);
-	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
-	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");

src/post-refinement.c

@@ -526,7 +526,7 @@ void pr_refine(struct image *image, const RefList *full)
 
 		max_shift = pr_iterate(image, full);
 
-		update_partialities(image, NULL, NULL, NULL);
+		update_partialities(image);
 
 		if ( verbose ) {
 			dev = guide_dev(image, full);

src/reflist.c

@@ -162,6 +162,34 @@ RefList *reflist_new()
 }
 
 
+/**
+ * reflection_new:
+ * @h: The h index of the new reflection
+ * @k: The k index of the new reflection
+ * @l: The l index of the new reflection
+ *
+ * Creates a new individual reflection.  You'll probably want to use
+ * add_refl_to_list() at some later point.
+ */
+Reflection *reflection_new(signed int h, signed int k, signed int l)
+{
+	return new_node(SERIAL(h, k, l));
+}
+
+
+/**
+ * reflection_free:
+ * @refl: The reflection to free.
+ *
+ * Destroys an individual reflection.
+ */
+void reflection_free(Reflection *refl)
+{
+	pthread_mutex_destroy(&refl->lock);
+	free(refl);
+}
+
+
 static void recursive_free(Reflection *refl)
 {
 	if ( refl->child[0] != NULL ) recursive_free(refl->child[0]);
@@ -169,8 +197,7 @@ static void recursive_free(Reflection *refl)
 
 	while ( refl != NULL ) {
 		Reflection *next = refl->next;
-		pthread_mutex_destroy(&refl->lock);
-		free(refl);
+		reflection_free(refl);
 		refl = next;
 	}
 }
@@ -806,6 +833,34 @@ Reflection *add_refl(RefList *list, signed int h, signed int k, signed int l)
 }
 
 
+/**
+ * add_refl_to_list
+ * @refl: A %Reflection
+ * @list: A %RefList
+ *
+ * Adds a reflection to @list.  The reflection that actually gets added will be
+ * a newly created one, and all the data will be copied across.  The original
+ * reflection will be destroyed and the new reflection returned.
+ *
+ * Returns: The newly created reflection, or NULL on failure.
+ *
+ **/
+Reflection *add_refl_to_list(Reflection *refl, RefList *list)
+{
+	signed int h, k, l;
+	Reflection *r_added;
+
+	get_indices(refl, &h, &k, &l);
+	r_added = add_refl(list, h, k, l);
+	if ( r_added == NULL ) return NULL;
+
+	copy_data(r_added, refl);
+	reflection_free(refl);
+
+	return r_added;
+}
+
+
 /********************************* Iteration **********************************/
 
 struct _reflistiterator {

src/reflist.h

@@ -50,6 +50,8 @@ typedef struct _reflistiterator RefListIterator;
 /* Creation/deletion */
 extern RefList *reflist_new(void);
 extern void reflist_free(RefList *list);
+extern Reflection *reflection_new(signed int h, signed int k, signed int l);
+extern void reflection_free(Reflection *refl);
 
 /* Search */
 extern Reflection *find_refl(const RefList *list, signed int h, signed int k, signed int l);
@@ -95,6 +97,7 @@ extern void set_symmetric_indices(Reflection *refl,
 /* Insertion */
 extern Reflection *add_refl(RefList *list,
                             signed int h, signed int k, signed int l);
+extern Reflection *add_refl_to_list(Reflection *refl, RefList *list);
 
 /* Iteration */
 extern Reflection *first_refl(RefList *list, RefListIterator **piter);