WIP on bringing programs up to date

8 files changed, 165 insertions, 75 deletions

diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c
index 3c0649cc..8ba955f6 100644
--- a/libcrystfel/src/crystal.c
+++ b/libcrystfel/src/crystal.c
@@ -55,7 +55,7 @@ struct _crystal
 	double                  osf;
 	double                  profile_radius;
 	int                     pr_dud;
-	double                  diffracting_resolution;
+	double                  resolution_limit;
 
 	/* Integrated (or about-to-be-integrated) reflections */
 	RefList                 *reflections;
@@ -63,7 +63,7 @@ struct _crystal
 };
 
 
-/************************** Setters and Constructors **************************/
+/************************** Constructor / Destructor **************************/
 
 
 /**
@@ -83,7 +83,7 @@ Crystal *crystal_new()
 
 	cryst->cell = NULL;
 	cryst->reflections = NULL;
-	cryst->diffracting_resolution = 0.0;
+	cryst->resolution_limit = 0.0;
 	cryst->n_saturated = 0;
 
 	return cryst;
@@ -97,11 +97,78 @@ Crystal *crystal_new()
  * Frees a %Crystal, and all internal resources concerning that crystal.
  *
  */
-void crystal_free(UnitCell *cryst)
+void crystal_free(Crystal *cryst)
 {
 	if ( cryst == NULL ) return;
-	free(crysta);
+	if ( cryst->cell != NULL ) cell_free(cryst->cell);
+	if ( cryst->reflections != NULL ) reflist_free(cryst->reflections);
+	free(cryst);
 }
 
 
 /********************************** Getters ***********************************/
+
+
+UnitCell *crystal_get_cell(Crystal *cryst)
+{
+	return cryst->cell;
+}
+
+
+double crystal_get_profile_radius(Crystal *cryst)
+{
+	return cryst->profile_radius;
+}
+
+
+RefList *crystal_get_reflections(Crystal *cryst)
+{
+	return cryst->reflections;
+}
+
+
+double crystal_get_resolution_limit(Crystal *cryst)
+{
+	return cryst->resolution_limit;
+}
+
+
+long long int crystal_get_num_saturated_reflections(Crystal *cryst)
+{
+	return cryst->n_saturated;
+}
+
+
+/********************************** Setters ***********************************/
+
+
+void crystal_set_cell(Crystal *cryst, UnitCell *cell)
+{
+	if ( cryst->cell != NULL ) cell_free(cryst->cell);
+	cryst->cell = cell_new_from_cell(cell);
+}
+
+
+void crystal_set_profile_radius(Crystal *cryst, double r)
+{
+	cryst->profile_radius = r;
+}
+
+
+void crystal_set_reflections(Crystal *cryst, RefList *reflist)
+{
+	if ( cryst->reflections != NULL ) reflist_free(reflist);
+	cryst->reflections = reflist;
+}
+
+
+void crystal_set_resolution_limit(Crystal *cryst, double res)
+{
+	cryst->resolution_limit = res;
+}
+
+
+void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n)
+{
+	cryst->n_saturated = n;
+}
diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h
index e5e40589..ca215399 100644
--- a/libcrystfel/src/crystal.h
+++ b/libcrystfel/src/crystal.h
@@ -60,7 +60,8 @@ extern void crystal_set_cell(Crystal *cryst, UnitCell *cell);
 extern void crystal_set_profile_radius(Crystal *cryst, double r);
 extern void crystal_set_reflections(Crystal *cryst, RefList *reflist);
 extern void crystal_set_resolution_limit(Crystal *cryst, double res);
-extern void crystal_set_num_saturated_reflections(Crystal *cryst, int n);
+extern void crystal_set_num_saturated_reflections(Crystal *cryst,
+                                                  long long int n);
 
 
 #endif	/* CRYSTAL_H */
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 0b2debda..6fbb9663 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -128,17 +128,20 @@ static int check_cell(struct dirax_private *dp, struct image *image,
 		return 0;
 	}
 
-	crystal_set_cell(cr, cell);
+	crystal_set_cell(cr, out);
 
 	if ( dp->indm & INDEXING_CHECK_PEAKS ) {
 		if ( !peak_sanity_check(image, &cr, 1) ) {
 			crystal_free(cr);  /* Frees the cell as well */
+			cell_free(out);
 			return 0;
 		}
 	}
 
 	image_add_crystal(image, cr);
 
+	cell_free(out); /* Crystal makes its own copy */
+
 	return 1;
 }
 
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 093e20f2..6e4fd40a 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -159,3 +159,21 @@ void image_remove_feature(ImageFeatureList *flist, int idx)
 {
 	flist->features[idx].valid = 0;
 }
+
+
+void image_add_crystal(struct image *image, Crystal *cryst)
+{
+	Crystal **crs;
+	int n;
+
+	n = image->n_crystals;
+	crs = realloc(image->crystals, (n+1)*sizeof(Crystal *));
+	if ( crs == NULL ) {
+		ERROR("Failed to allocate memory for crystals.\n");
+		return;
+	}
+
+	crs[n] = cryst;
+	image->crystals = crs;
+	image->n_crystals = n+1;
+}
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 8b300012..d1d3c134 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -149,7 +149,7 @@ static int check_cell(struct mosflm_private *mp, struct image *image,
 		return 0;
 	}
 
-	crystal_set_cell(cr, cell);
+	crystal_set_cell(cr, out);
 
 	if ( mp->indm & INDEXING_CHECK_PEAKS ) {
 		if ( !peak_sanity_check(image, &cr, 1) ) {
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 7529d12f..b144e266 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -103,7 +103,7 @@ struct reax_search
 
 struct reax_private
 {
-	IndexingPrivate base;
+	IndexingMethod indm;
 	struct dvec *directions;
 	int n_dir;
 	double angular_inc;
@@ -915,9 +915,39 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew,
 }
 
 
+static int check_cell(struct reax_private *dp, struct image *image,
+                      UnitCell *cell)
+{
+	UnitCell *out;
+	Crystal *cr;
+
+	out = cell_new_from_cell(cell);
+
+	cr = crystal_new();
+	if ( cr == NULL ) {
+		ERROR("Failed to allocate crystal.\n");
+		return 0;
+	}
+
+	crystal_set_cell(cr, out);
+
+	if ( dp->indm & INDEXING_CHECK_PEAKS ) {
+		if ( !peak_sanity_check(image, &cr, 1) ) {
+			crystal_free(cr);  /* Frees the cell as well */
+			return 0;
+		}
+	}
+
+	image_add_crystal(image, cr);
+
+	return 1;
+}
+
+
 static void assemble_cells_from_candidates(struct image *image,
                                            struct reax_search *s,
-                                           UnitCell *cell)
+                                           UnitCell *cell,
+                                           struct reax_private *p)
 {
 	int i, j, k;
 	signed int ti, tj, tk;
@@ -967,7 +997,9 @@ static void assemble_cells_from_candidates(struct image *image,
 			continue;
 		}
 
-		peak_lattice_agreement(image, cnew, &fom);
+		/* FIXME! */
+		//peak_lattice_agreement(image, cnew, &fom);
+		fom = 1.0;
 		add_cell_candidate(&cl, cnew, fom);
 
 	}
@@ -985,22 +1017,20 @@ static void assemble_cells_from_candidates(struct image *image,
 		cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga);
 		cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA,
 		                                     &alA, &beA, &gaA);
-		if ( (a - aA) > aA/10.0 ) w = 1;
-		if ( (b - bA) > bA/10.0 ) w = 1;
-		if ( (c - cA) > cA/10.0 ) w = 1;
-		if ( (al - alA) > deg2rad(5.0) ) w = 1;
-		if ( (be - beA) > deg2rad(5.0) ) w = 1;
-		if ( (ga - gaA) > deg2rad(5.0) ) w = 1;
+		if ( fabs(a - aA) > aA/10.0 ) w = 1;
+		if ( fabs(b - bA) > bA/10.0 ) w = 1;
+		if ( fabs(c - cA) > cA/10.0 ) w = 1;
+		if ( fabs(al - alA) > deg2rad(5.0) ) w = 1;
+		if ( fabs(be - beA) > deg2rad(5.0) ) w = 1;
+		if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1;
 		if ( w ) {
 			STATUS("This cell is a long way from that sought:\n");
 			cell_print(cl.cand[i].cell);
 		}
 	}
 
-	image->ncells = cl.n_cand;
-	assert(image->ncells <= MAX_CELL_CANDIDATES);
 	for ( i=0; i<cl.n_cand; i++ ) {
-		image->candidate_cells[i] = cl.cand[i].cell;
+		if ( check_cell(p,  image, cl.cand[i].cell) ) break;
 	}
 
 	free(cl.cand);
@@ -1016,7 +1046,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 	struct reax_search *s;
 	int i;
 
-	assert(pp->indm == INDEXING_REAX);
 	p = (struct reax_private *)pp;
 
 	fft_in = fftw_malloc(p->nel*sizeof(double));
@@ -1039,7 +1068,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 //	                        fft_in, fft_out, p->plan, smin, smax,
 //	                        image->det, p);
 
-	assemble_cells_from_candidates(image, s, cell);
+	assemble_cells_from_candidates(image, s, cell, p);
 
 	for ( i=0; i<s->n_search; i++ ) {
 		free(s->search[i].cand);
@@ -1051,7 +1080,9 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 }
 
 
-IndexingPrivate *reax_prepare()
+IndexingPrivate *reax_prepare(IndexingMethod indm, UnitCell *cell,
+                              const char *filename, struct detector *det,
+                              struct beam_params *beam, float *ltl)
 {
 	struct reax_private *p;
 	int samp;
@@ -1060,8 +1091,6 @@ IndexingPrivate *reax_prepare()
 	p = calloc(1, sizeof(*p));
 	if ( p == NULL ) return NULL;
 
-	p->base.indm = INDEXING_REAX;
-
 	p->angular_inc = deg2rad(1.0);
 
 	/* Reserve memory, over-estimating the number of directions */
@@ -1122,6 +1151,8 @@ IndexingPrivate *reax_prepare()
 	p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out,
 	                             1, FFTW_MEASURE);
 
+	p->indm = indm;
+
 	return (IndexingPrivate *)p;
 }
 
@@ -1130,7 +1161,6 @@ void reax_cleanup(IndexingPrivate *pp)
 {
 	struct reax_private *p;
 
-	assert(pp->indm == INDEXING_REAX);
 	p = (struct reax_private *)pp;
 
 	free(p->directions);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 943ee656..31fb8527 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -69,30 +69,6 @@ struct _stream
 };
 
 
-int count_patterns(FILE *fh)
-{
-	char *rval;
-
-	int n_total_patterns = 0;
-	do {
-		char line[1024];
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++;
-
-	} while ( rval != NULL );
-
-	if ( ferror(fh) ) {
-		ERROR("Read error while counting patterns.\n");
-		return 0;
-	}
-
-	return n_total_patterns;
-}
-
-
 static int read_peaks(FILE *fh, struct image *image)
 {
 	char *rval = NULL;
@@ -482,28 +458,6 @@ void write_stream_header(FILE *ofh, int argc, char *argv[])
 }
 
 
-int skip_some_files(FILE *fh, int n)
-{
-	char *rval = NULL;
-	int n_patterns = 0;
-
-	do {
-
-		char line[1024];
-
-		if ( n_patterns == n ) return 0;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
-
-	} while ( rval != NULL );
-
-	return 1;
-}
-
-
 Stream *open_stream_for_read(const char *filename)
 {
 	Stream *st;
@@ -591,3 +545,21 @@ int is_stream(const char *filename)
 
 	return 0;
 }
+
+
+/**
+ * rewind_stream:
+ * @st: A %Stream
+ *
+ * Attempts to set the file pointer for @st to the start of the stream, so that
+ * later calls to read_chunk() will repeat the sequence of chunks from the
+ * start.
+ *
+ * Programs must not assume that this operation always succeeds!
+ *
+ * Returns: non-zero if the stream could not be rewound.
+ */
+int rewind_stream(Stream *st)
+{
+	return fseek(st->fh, 0, SEEK_SET);
+}
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index db74583d..9340ca0d 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -49,9 +49,8 @@ extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile,
 
 extern int read_chunk(Stream *st, struct image *image);
 
-/* Nasty functions that should be avoided */
-extern int count_patterns(FILE *fh);
-extern int skip_some_files(FILE *fh, int n);
+extern int rewind_stream(Stream *st);
+
 extern int is_stream(const char *filename);
 
 #endif	/* STREAM_H */