Removed global variable, include profile_radius and diff_limit in sol file, cleaned up warning messages, not working yet for multi-crystals

1 files changed, 145 insertions, 120 deletions

diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 2b5907ba..0150787a 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -42,6 +42,7 @@
 #include <string.h>
 #include <assert.h>
 #include <fenv.h>
+#include <unistd.h>
 
 #include "image.h"
 #include "utils.h"
@@ -65,28 +66,28 @@
 
 /** \file index.h */
 
+#define NPARAMS_PER_LINE 13  // There are 9 vector components, 2 detector shifts, 1 profile radius, 1 resolution limit 
+#define NKEYS_PER_LINE 4  // The keys are the filename, event path, event dim and crystal number 
+
 struct skip_private
 {
-	char path_to_sol[50];
-    UnitCell *cellTemplate;
-	float solutions[];
+	UnitCell *cellTemplate;
+	struct record_t *sol_hash;
 };
 
-struct record_key_t{			 	//added
+struct record_key_t{			 	//Hash table keys
   char filename[100];
-  char event_path[35];
+  char event_path[100];
   int event_dim;
+  int crystal_number;
 };
 
-struct record_t{			 	//added
+struct record_t{			 		//Hash table
     struct record_key_t key;
-    int line;
+    float solution[NPARAMS_PER_LINE];
     UT_hash_handle hh;
 };
 
-//declared as global variable
-struct record_t *line_hash = NULL;		//added
-
 struct _indexingprivate
 {
 	IndexingFlags flags;
@@ -98,22 +99,24 @@ struct _indexingprivate
 	void **engine_private;
 };
 
-void delete_all() {
-  struct record_t *current_line, *tmp;
+void print_struct(struct record_t *sol_hash) {
+    struct record_t *s;
+    s = (struct record_t *)malloc(sizeof *s);
+    memset(s, 0, sizeof *s);
 
-  HASH_ITER(hh, line_hash, current_line, tmp) {
-    HASH_DEL(line_hash, current_line);
-    free(current_line);             
-  }
+    for(s=sol_hash; s != NULL; s=(struct record_t*)(s->hh.next)) {
+        printf("This entry corresponds to event %s//%d and crystal_number %d in file %s \n", s->key.event_path, s->key.event_dim, s->key.crystal_number, s->key.filename);
+    }
 }
 
-void print_struct() {
+void full_print_struct(struct record_t *sol_hash) {
     struct record_t *s;
     s = (struct record_t *)malloc(sizeof *s);
     memset(s, 0, sizeof *s);
 
-    for(s=line_hash; s != NULL; s=(struct record_t*)(s->hh.next)) {
-        printf("The line is event_path %s and event_dim %d in file %s \n", s->key.event_path, s->key.event_dim, s->key.filename);
+    for(s=sol_hash; s != NULL; s=(struct record_t*)(s->hh.next)) {
+        printf("This entry corresponds to event %s//%d and crystal_number %d in file %s \n", s->key.event_path, s->key.event_dim, s->key.crystal_number, s->key.filename);
+		printf("The solution parameters are %e %e %e %e %e %e %e %e %e %e %e %e %e \n", s->solution[0], s->solution[1], s->solution[2], s->solution[3], s->solution[4], s->solution[5], s->solution[6], s->solution[7], s->solution[8], s->solution[9], s->solution[10], s->solution[11], s->solution[12]);
     }
 }
 
@@ -137,7 +140,6 @@ static void set_last_task(char *lt, const char *task)
 	strcpy(lt, task);
 }
 
-
 static void show_indexing_flags(IndexingFlags flags)
 {
 	STATUS("Indexing parameters:\n");
@@ -153,30 +155,60 @@ static void show_indexing_flags(IndexingFlags flags)
 	       onoff(flags & INDEXING_RETRY));
 }
 
+int ncrystals_in_sol(char *path)
+{
+	FILE *fh;
+	int count = 0;  // Line counter (result) 
+    char c;  // To store a character read from file 
+
+	fh = fopen(path, "r");
+
+	if ( fh == NULL ) {
+		ERROR("%s not found by ncrystals_in_sol\n",path);
+		return 0;
+	}
+
+	for (c = getc(fh); c != EOF; c = getc(fh)) 
+        if (c == '\n') // Increment count if this character is newline 
+            count = count + 1; 
+	
+	//For the last line, which has no \n at the end
+	count = count + 1;
+  
+    // Close the file 
+    fclose(fh); 
+
+	// STATUS("Found indexing file %s containing %d lines. \n", path, count);
+
+	return count;
+
+}
+
 
 void *skip_prepare(char *solution_filename, UnitCell *cell)
 {	
 	FILE *fh;
 	int nlines;
 	int nparams_in_solution;
-	int nparams_per_line;
 	int nentries;
-	char *filename[35];   					//??? is that correct
-	char *event_path[35];                   //??? is that correct
-	int event_dim;
+	char filename[100];   					
+	char event_path[100];   
+	int event_dim;                
+	int crystal_number;
 	int current_line;
 	int position_in_current_line;
-	struct record_t *item;
+	struct record_t *sol_hash = NULL;
+	struct record_t *item = NULL;
+	float params[NPARAMS_PER_LINE];
+	char path_to_sol[50],extension[10];
+	char cwd[PATH_MAX];
 
 	//Assembling solution file name from input file name to crystfel
-	char path_to_sol[50],extension[10];
 	strcpy(path_to_sol,"../");
 	strcpy(extension,".sol");
-	strtok(solution_filename, "."); //Takes the first part of the list name (before the .)
 	strcat(path_to_sol, strtok(solution_filename, ".") );     //Add .sol at the previously splitted name
 	strcat(path_to_sol, extension );
 
-	char cwd[PATH_MAX];
 	if (getcwd(cwd, sizeof(cwd)) != NULL) {
 		//printf("Path where skip_prepare is ran: %s\n", cwd);
 	}
@@ -192,22 +224,19 @@ void *skip_prepare(char *solution_filename, UnitCell *cell)
 	}
 
 	nlines = ncrystals_in_sol(path_to_sol);
-	nparams_per_line = 11;  /* There are 9 vector components and 2 detector shifts */
-	nparams_in_solution = nlines*nparams_per_line; /* total crystal parameters in solution file */
-	nentries = nlines*(nparams_per_line+3); /* total entries in solution file */
+	nparams_in_solution = nlines*NPARAMS_PER_LINE; /* total crystal parameters in solution file */
+	nentries = nlines*(NPARAMS_PER_LINE+NKEYS_PER_LINE); /* total entries in solution file */
   
-	float *params = malloc(nparams_in_solution * sizeof( float));
-
 	STATUS("Parsing solution file containing %d lines...\n", nlines);
 
 	//Reads indexing solutions
-	int j = 0; //index that follows the solution parameter
+	int j = 0; //index that follows the solution parameter [0,NPARAMS_PER_LINE]
 	for(int i = 0; i < nentries; i++)
 	{	
 
-		current_line = i/(nparams_per_line+3);
+		current_line = i/(NPARAMS_PER_LINE+NKEYS_PER_LINE);
 		
-		position_in_current_line = (i)%(nparams_per_line+3);
+		position_in_current_line = (i)%(NPARAMS_PER_LINE+NKEYS_PER_LINE);
 
 		if (position_in_current_line == 0){
 
@@ -225,106 +254,84 @@ void *skip_prepare(char *solution_filename, UnitCell *cell)
 		if (position_in_current_line == 1){
 
 			if ( fscanf(fh, "%s", event_path) != 1 ) {
-				printf("Failed to read an event_path from indexing.sol with skip_prepare\n");
+				printf("Failed to read an event path from indexing.sol with skip_prepare\n");
 				return 0;
 			}
 		}
-			
+
 		if (position_in_current_line == 2){
+
 			if ( fscanf(fh, "%d", &event_dim) != 1 ) {
-				printf("Failed to read an event from indexing.sol with skip_prepare\n");
+				printf("Failed to read an event dim from indexing.sol with skip_prepare\n");
+				return 0;
+			}
+		}
+			
+		if (position_in_current_line == 3){
+			if ( fscanf(fh, "%d", &crystal_number) != 1 ) {
+				printf("Failed to read a crystal number from indexing.sol with skip_prepare\n");
+				return 0;
+			}			
+		}
+
+		if (position_in_current_line > 3){
+			if ( fscanf(fh, "%e", &params[j]) != 1 ) {
+				printf("Failed to read a parameter from indexing.sol with skip_prepare\n");
 				return 0;
 			}
+			j+=1;
+
+			
+		}
+
+		if (j == (NPARAMS_PER_LINE) ){
 
-			//Prepare to add to the hash table once an event is read
+			//Prepare to add to the hash table
 			item = (struct record_t *)malloc(sizeof *item);
 			memset(item, 0, sizeof *item);
 			strcpy(item->key.filename, filename);
 			strcpy(item->key.event_path, event_path);
 			item->key.event_dim = event_dim;
-			item->line=current_line;
+			item->key.crystal_number = crystal_number;
+			for (int k = 0; k < NPARAMS_PER_LINE; k++){
+    			item->solution[k] = params[k];
+			}
 
 			//Verify the uniqueness of the key
 			struct record_t *uniqueness_test;
-			HASH_FIND(hh, line_hash, &item->key, sizeof(struct record_key_t), uniqueness_test);  /* id already in the hash? */
+			HASH_FIND(hh, sol_hash, &item->key, sizeof(struct record_key_t), uniqueness_test);  // id already in the hash? 
     		if (uniqueness_test==NULL) {
 				//If not in the table, add
-				HASH_ADD(hh, line_hash, key, sizeof(struct record_key_t), item);
+				HASH_ADD(hh, sol_hash, key, sizeof(struct record_key_t), item);
 			}
 			else{
 				//If already in the table, break with error
-				printf("Keys to the data must be unique! Verify the combination filename, event_path, event_dim");
+				printf("Keys to the data must be unique! Verify the combination filename, event, crystal_number");
 				return 0;
 			}
 
+		j=0;		
 
-			
-		}
-
-		if (position_in_current_line > 2){
-			if ( fscanf(fh, "%e", &params[j]) != 1 ) {
-				printf("Failed to read a parameter from indexing.sol with skip_prepare\n");
-				return 0;
-			}
-			j+=1;
 		}
 	}
 	
 	fclose(fh);
-
+	
 	STATUS("Solution file parsing done. Have %d parameters and %d total entries.\n", nparams_in_solution, nentries)
 
 	struct skip_private *dp;
-	dp = (struct skip_private *) malloc( sizeof(struct skip_private) + nentries * sizeof( float) );
+	dp = (struct skip_private *) malloc( sizeof(struct skip_private));
 
 	if ( dp == NULL ) return NULL;
-	
-	memcpy(dp->path_to_sol, path_to_sol, sizeof path_to_sol);
-
-	for (int k = 0; k < nparams_in_solution; k++)
-	{
-		dp->solutions[k] = params[k];
-	}
-    
     
     dp->cellTemplate = cell;
-
-	free(params);
+	dp->sol_hash = sol_hash;
 	
 	STATUS("Solution lookup table initialized!\n");
 
 	return (void *)dp;
 }
 
-int ncrystals_in_sol(char *path)
-{
-	FILE *fh;
-	int count = 0;  // Line counter (result) 
-    char c;  // To store a character read from file 
-
-	fh = fopen(path, "r");
-
-	if ( fh == NULL ) {
-		ERROR("%s not found by ncrystals_in_sol\n",path);
-		return 0;
-	}
-
-	for (c = getc(fh); c != EOF; c = getc(fh)) 
-        if (c == '\n') // Increment count if this character is newline 
-            count = count + 1; 
-	
-	//For the last line, which has no \n at the end
-	count = count + 1;
-  
-    // Close the file 
-    fclose(fh); 
-
-	// STATUS("Found indexing file %s containing %d lines. \n", path, count);
-
-	return count;
-
-}
-
 //Fonction from pinkindexer.c to update the detector center of individual crystals
 //hack for electron crystallography while crystal_set_det_shift is not working approprietly
 static void update_detector(struct detector *det, double xoffs, double yoffs)
@@ -338,15 +345,13 @@ static void update_detector(struct detector *det, double xoffs, double yoffs)
 	}
 }
 
-static int skip_index(struct image *image, void *mpriv)
+static int skip_index(struct image *image, void *mpriv, int crystal_number)
 {
 	Crystal *cr;
 	UnitCell *cell;
-	float asx, asy, asz, bsx, bsy, bsz, csx, csy, csz, xshift, yshift;
-	int serial;
-	int nparams_per_line;
-	struct record_t *item, *p;
-	int line_number;
+	float asx, asy, asz, bsx, bsy, bsz, csx, csy, csz, xshift, yshift, profile_radius, resolution_limit;
+	struct record_t *item, *p, *pprime;
+	float *sol;
 	
 	struct skip_private *dp = (struct skip_private *)mpriv;
 
@@ -356,29 +361,29 @@ static int skip_index(struct image *image, void *mpriv)
 	strcpy(item->key.filename, image->filename);
 	strcpy(item->key.event_path, *image->event->path_entries);
 	item->key.event_dim = *image->event->dim_entries;
+	item->key.crystal_number = crystal_number;
 
-    HASH_FIND(hh, line_hash, &item->key, sizeof(struct record_key_t), p);  /* id already in the hash? */
+    HASH_FIND(hh,  dp->sol_hash, &item->key, sizeof(struct record_key_t), p);  //id already in the hash?
     if (p==NULL) {
-		/*
-		STATUS("Not indexing file %s, event %s %d \n", image->filename, *image->event->path_entries, *image->event->dim_entries);
-		*/
+		//STATUS("Not indexing file %s, event %s %d \n", image->filename, *image->event->path_entries, *image->event->dim_entries);
 		return 0;
 	}
 
-	line_number = p->line;
-	nparams_per_line = 11;
-
-	asx = dp->solutions[(line_number * nparams_per_line) + 0];
-	asy = dp->solutions[(line_number * nparams_per_line) + 1];
-	asz = dp->solutions[(line_number * nparams_per_line) + 2];
-	bsx = dp->solutions[(line_number * nparams_per_line) + 3];
-	bsy = dp->solutions[(line_number * nparams_per_line) + 4];
-	bsz = dp->solutions[(line_number * nparams_per_line) + 5];
-	csx = dp->solutions[(line_number * nparams_per_line) + 6];
-	csy = dp->solutions[(line_number * nparams_per_line) + 7];
-	csz = dp->solutions[(line_number * nparams_per_line) + 8];
-	xshift = dp->solutions[(line_number * nparams_per_line) + 9];
-	yshift = dp->solutions[(line_number * nparams_per_line) + 10];
+	sol = &(p->solution)[0];
+	
+	asx = sol[0];
+	asy = sol[1];
+	asz = sol[2];
+	bsx = sol[3];
+	bsy = sol[4];
+	bsz = sol[5];
+	csx = sol[6];
+	csy = sol[7];
+	csz = sol[8];
+	xshift = sol[9];
+	yshift = sol[10];
+	profile_radius = sol[11];
+	resolution_limit = sol[12];
 
 	cr = crystal_new();
 	cell = cell_new();
@@ -390,9 +395,29 @@ static int skip_index(struct image *image, void *mpriv)
 	crystal_set_cell(cr, cell);
 	crystal_set_det_shift(cr, xshift * 1e-3, yshift * 1e-3);
 	update_detector(image->det, xshift * 1e-3, yshift * 1e-3);
+	crystal_set_profile_radius(cr, profile_radius*1e9);
+	crystal_set_resolution_limit(cr, resolution_limit*1e9);
 	image_add_crystal(image, cr);
 
+	//Look for additional crystals
+	item->key.crystal_number = crystal_number+1;
+	HASH_FIND(hh,  dp->sol_hash, &item->key, sizeof(struct record_key_t), pprime);  //id already in the hash?
+    if (pprime==NULL) {
+		
+		//If no more crystal, done
+		return 1;
+	}
+	else{
+		//If more crystals, recursive call for next crystal in line
+		skip_index(image,mpriv,crystal_number+1);
+	}
+
+	
+	dp->sol_hash = NULL; //Clean up local copy of the hash table
+	
+
 	return 1;
+	
 }
 
 
@@ -771,7 +796,6 @@ void cleanup_indexing(IndexingPrivate *ipriv)
 
 			case INDEXING_FILE :
 			free(ipriv->engine_private[n]);
-			delete_all();
 			break;
 
 			case INDEXING_TAKETWO :
@@ -899,7 +923,8 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 
 		case INDEXING_FILE :
 		set_last_task(last_task, "indexing:file");
-		r = skip_index(image,mpriv);
+		int crystal_number = 0;
+		r = skip_index(image,mpriv,crystal_number);
 		break;
 
 		case INDEXING_FELIX :