Put ReAx cell through reduction, and add initial panel refinement stuff

1 files changed, 114 insertions, 11 deletions

diff --git a/src/reax.c b/src/reax.c
index eafb1225..0af735ac 100644
--- a/src/reax.c
+++ b/src/reax.c
@@ -61,7 +61,8 @@ struct reax_private
 static double check_dir(struct dvec *dir, ImageFeatureList *flist,
                         int nel, double pmax, double *fft_in,
                         fftw_complex *fft_out, fftw_plan plan,
-                        int smin, int smax)
+                        int smin, int smax,
+                        const char *rg, struct detector *det)
 {
 	int n, i;
 	double tot;
@@ -80,6 +81,17 @@ static double check_dir(struct dvec *dir, ImageFeatureList *flist,
 		f = image_get_feature(flist, i);
 		if ( f == NULL ) continue;
 
+		if ( rg != NULL ) {
+
+			struct panel *p;
+
+			p = find_panel(det, f->fs, f->ss);
+			assert(p != NULL);
+
+			if ( p->rigid_group != rg ) continue;
+
+		}
+
 		val = f->rx*dir->x + f->ry*dir->y + f->rz*dir->z;
 
 		idx = nel/2 + nel*val/(2.0*pmax);
@@ -245,7 +257,7 @@ static void fine_search(struct reax_private *p, ImageFeatureList *flist,
 		dir.z = cos(th);
 
 		new_fom = check_dir(&dir, flist, nel, pmax, fft_in, fft_out,
-		                    plan, smin, smax);
+		                    plan, smin, smax, NULL, NULL);
 
 		if ( new_fom > fom ) {
 			fom = new_fom;
@@ -277,6 +289,93 @@ static void fine_search(struct reax_private *p, ImageFeatureList *flist,
 }
 
 
+static double get_model_phase(double x, double y, double z, ImageFeatureList *f,
+                              int nel, double pmax, double *fft_in,
+                              fftw_complex *fft_out, fftw_plan plan,
+                              int smin, int smax, const char *rg,
+                              struct detector *det)
+{
+	struct dvec dir;
+	int s, i;
+	double max;
+	double re, im;
+
+	dir.x = x;  dir.y = y;  dir.z = z;
+
+	check_dir(&dir, f, nel, pmax, fft_in,fft_out, plan, smin, smax,
+	          rg, det);
+
+	s = -1;
+	max = 0.0;
+	for ( i=smin; i<=smax; i++ ) {
+
+		double re, im, m;
+
+		re = fft_out[i][0];
+		im = fft_out[i][1];
+		m = sqrt(re*re + im*im);
+		if ( m > max ) {
+			max = m;
+			s = i;
+		}
+
+	}
+
+	re = fft_out[s][0];
+	im = fft_out[s][1];
+
+	return atan2(im, re);
+}
+
+
+static void refine_rigid_group(struct image *image, UnitCell *cell,
+                               const char *rg, int nel, double pmax,
+                               double *fft_in, fftw_complex *fft_out,
+                               fftw_plan plan, int smin, int smax,
+                               struct detector *det)
+{
+	double ax, ay, az, ma;
+	double bx, by, bz, mb;
+	double cx, cy, cz, mc;
+	double pha, phb, phc;
+
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	ma = modulus(ax, ay, az);
+	mb = modulus(bx, by, bz);
+	mc = modulus(cx, cy, cz);
+
+	ax /= ma;  ay /= ma;  az /= ma;
+	bx /= mb;  by /= mb;  bz /= mb;
+	cx /= mc;  cy /= mc;  cz /= mc;
+
+	pha = get_model_phase(ax, ay, az, image->features, nel, pmax, fft_in,
+	                      fft_out, plan, smin, smax, rg, det);
+	phb = get_model_phase(bx, by, bz, image->features, nel, pmax, fft_in,
+	                      fft_out, plan, smin, smax, rg, det);
+	phc = get_model_phase(cx, cy, cz, image->features, nel, pmax, fft_in,
+	                      fft_out, plan, smin, smax, rg, det);
+
+
+}
+
+
+static void refine_all_rigid_groups(struct image *image, UnitCell *cell,
+                                    int nel, double pmax,
+                                    double *fft_in, fftw_complex *fft_out,
+                                    fftw_plan plan, int smin, int smax,
+                                    struct detector *det)
+{
+	int i;
+
+	for ( i=0; i<image->det->num_rigid_groups; i++ ) {
+		refine_rigid_group(image, cell, image->det->rigid_groups[i],
+		                   nel, pmax, fft_in, fft_out, plan, smin, smax,
+		                   det);
+	}
+}
+
+
 void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 {
 	int i;
@@ -352,7 +451,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 
 		new_fom = check_dir(&p->directions[i], image->features,
 		                    p->nel, pmax, fft_in, fft_out, p->plan,
-		                    smin, smax);
+		                    smin, smax, NULL, NULL);
 		if ( new_fom > fom ) {
 			fom = new_fom;
 			th = p->directions[i].th;
@@ -377,7 +476,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 
 		new_fom = check_dir(&p->directions[i], image->features,
 		                    p->nel, pmax, fft_in, fft_out, p->plan,
-		                    smin, smax);
+		                    smin, smax, NULL, NULL);
 		if ( new_fom > fom ) {
 			fom = new_fom;
 			th = p->directions[i].th;
@@ -406,7 +505,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 
 		new_fom = check_dir(&p->directions[i], image->features,
 		                    p->nel, pmax, fft_in, fft_out, p->plan,
-		                    smin, smax);
+		                    smin, smax, NULL, NULL);
 		if ( new_fom > fom ) {
 			fom = new_fom;
 			th = p->directions[i].th;
@@ -417,16 +516,20 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell)
 	fine_search(p, image->features, p->nel, pmax, fft_in, fft_out,
 	            p->plan, smin, smax, th, ph, &cx, &cy, &cz);
 
-	image->indexed_cell = cell_new();
-	cell_set_cartesian(image->indexed_cell, ax, ay, az,
-	                                        bx, by, bz, cx, cy, cz);
-	refine_cell(image, image->indexed_cell, image->features);
+	image->candidate_cells[0] = cell_new();
+	cell_set_cartesian(image->candidate_cells[0],
+	                   ax, ay, az, bx, by, bz, cx, cy, cz);
+
+	refine_all_rigid_groups(image, image->candidate_cells[0], p->nel, pmax,
+	                        fft_in, fft_out, p->plan, smin, smax,
+	                        image->det);
+	map_all_peaks(image);
+	refine_cell(image, image->candidate_cells[0], image->features);
 
 	fftw_free(fft_in);
 	fftw_free(fft_out);
 
-	/* We fill in indexed_cell but none of the candidate_cells, which
-	 * bypasses any cell reduction which is selected. */
+	image->ncells = 1;
 }