aboutsummaryrefslogtreecommitdiff
path: root/src/whirligig.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/whirligig.c')
-rw-r--r--src/whirligig.c144
1 files changed, 119 insertions, 25 deletions
diff --git a/src/whirligig.c b/src/whirligig.c
index a676d707..ed25e887 100644
--- a/src/whirligig.c
+++ b/src/whirligig.c
@@ -44,28 +44,114 @@
#include <stream.h>
#include "version.h"
+#include "cell-utils.h"
+#include "integer_matrix.h"
-static void process_series(struct image *images, signed int *ser, int len)
+static void process_series(struct image *images, signed int *ser,
+ IntegerMatrix **mat, int len)
{
- STATUS("Found a rotation series of %i views\n", len);
+// STATUS("Found a rotation series of %i views\n", len);
}
-static int gatinator(UnitCell *a, UnitCell *b)
+static int cells_are_similar(UnitCell *cell1, UnitCell *cell2)
{
+ double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1;
+ double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2;
+ UnitCellTransformation *tfn1, *tfn2;
+ UnitCell *pcell1, *pcell2;
+ const double atl = deg2rad(5.0);
+
+ /* Compare primitive cells, not centered */
+ pcell1 = uncenter_cell(cell1, &tfn1);
+ pcell2 = uncenter_cell(cell2, &tfn2);
+
+ cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1,
+ &bsx1, &bsy1, &bsz1,
+ &csx1, &csy1, &csz1);
+
+ cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2,
+ &bsx2, &bsy2, &bsz2,
+ &csx2, &csy2, &csz2);
+
+
+ cell_free(pcell1);
+ cell_free(pcell2);
+
+ if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0;
+ if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0;
+ if ( angle_between(csx1, csy1, csz1, csx2, csy2, csz2) > atl ) return 0;
+
+ /* FIXME: Moduli need to be similar as well */
+
+ return 1;
+}
+
+
+static int gatinator(UnitCell *a, UnitCell *b, IntegerMatrix **pmb)
+{
+ IntegerMatrix *m;
+ int i[9];
+
+ m = intmat_new(3, 3);
+
+ for ( i[0]=-1; i[0]<=+1; i[0]++ ) {
+ for ( i[1]=-1; i[1]<=+1; i[1]++ ) {
+ for ( i[2]=-1; i[2]<=+1; i[2]++ ) {
+ for ( i[3]=-1; i[3]<=+1; i[3]++ ) {
+ for ( i[4]=-1; i[4]<=+1; i[4]++ ) {
+ for ( i[5]=-1; i[5]<=+1; i[5]++ ) {
+ for ( i[6]=-1; i[6]<=+1; i[6]++ ) {
+ for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
+ for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
+
+ UnitCellTransformation *tfn;
+ UnitCell *nc;
+ int j, k;
+ int l = 0;
+
+ for ( j=0; j<3; j++ )
+ for ( k=0; k<3; k++ )
+ intmat_set(m, j, k, i[l++]);
+
+ if ( intmat_det(m) != +1 ) continue;
+
+ tfn = tfn_from_intmat(m);
+ nc = cell_transform(b, tfn);
+
+ if ( cells_are_similar(a, nc) ) {
+ *pmb = m;
+ return 1;
+ }
+
+ tfn_free(tfn);
+ cell_free(nc);
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ intmat_free(m);
return 0;
}
-static int try_all(struct image *a, struct image *b, int *c1, int *c2)
+static int try_all(struct image *a, struct image *b, int *c1, int *c2,
+ IntegerMatrix **m2)
{
int i, j;
for ( i=0; i<a->n_crystals; i++ ) {
for ( j=0; j<b->n_crystals; j++ ) {
if ( gatinator(crystal_get_cell(a->crystals[i]),
- crystal_get_cell(b->crystals[j])) )
+ crystal_get_cell(b->crystals[j]), m2) )
{
*c1 = i;
*c2 = j;
@@ -78,16 +164,19 @@ static int try_all(struct image *a, struct image *b, int *c1, int *c2)
}
-static void dump(struct image *win, signed int *series, int window_len, int pos)
+static void dump(struct image *win, signed int *ser, IntegerMatrix **mat,
+ int window_len, int pos)
{
int i;
for ( i=0; i<pos; i++ ) {
free_all_crystals(&win[i]);
+ intmat_free(mat[i]);
}
memmove(win, &win[pos], (window_len-pos)*sizeof(struct image *));
- memmove(series, &series[pos], (window_len-pos)*sizeof(signed int));
+ memmove(ser, &ser[pos], (window_len-pos)*sizeof(signed int));
+ memmove(mat, &mat[pos], (window_len-pos)*sizeof(IntegerMatrix *));
}
@@ -110,7 +199,8 @@ int main(int argc, char *argv[])
int polarisation = 1;
int pos = 0;
struct image *win;
- signed int *series;
+ signed int *ser;
+ IntegerMatrix **mat;
int window_len = 64;
/* Long options */
@@ -164,8 +254,9 @@ int main(int argc, char *argv[])
}
win = calloc(window_len, sizeof(struct image));
- series = calloc(window_len, sizeof(int));
- if ( (win == NULL) || (series == NULL) ) {
+ ser = calloc(window_len, sizeof(int));
+ mat = calloc(window_len, sizeof(IntegerMatrix *));
+ if ( (win == NULL) || (ser == NULL) || (mat == NULL) ) {
ERROR("Failed to allocate series buffers\n");
return 1;
}
@@ -193,30 +284,31 @@ int main(int argc, char *argv[])
/* Need at least two images to compare */
if ( pos == 0 ) {
- series[0] = -1;
+ ser[0] = -1;
pos++;
continue;
}
- if ( try_all(&win[pos-1], cur, &c1, &c2) ) {
- series[pos-1] = c1;
- series[pos] = c2;
- printf("-");
+ if ( try_all(&win[pos-1], cur, &c1, &c2, &mat[pos]) ) {
+ ser[pos-1] = c1;
+ ser[pos] = c2;
+ printf("*");
} else {
- series[pos] = -1;
- printf(".");
+ ser[pos] = -1;
+ mat[pos] = NULL;
+ printf("-");
}
fflush(stdout);
- if ( series[0] == -1 ) {
- dump(win, series, window_len, 1);
+ if ( ser[0] == -1 ) {
+ dump(win, ser, mat, window_len, 1);
pos--;
}
- if ( (series[pos] == -1) && (series[pos-1] == -1) ) {
+ if ( (ser[pos] == -1) && (ser[pos-1] == -1) ) {
/* Series ready to process */
- process_series(win, series, pos-2);
- dump(win, series, window_len, pos);
+ process_series(win, ser, mat, pos-2);
+ dump(win, ser, mat, window_len, pos);
pos = 0;
}
@@ -224,8 +316,9 @@ int main(int argc, char *argv[])
if ( pos == window_len ) {
window_len *= 2;
win = realloc(win, window_len*sizeof(struct image));
- series = realloc(series, window_len*sizeof(signed int));
- if ( (win == NULL) || (series == NULL) ) {
+ ser = realloc(ser, window_len*sizeof(signed int));
+ mat = realloc(mat, window_len*sizeof(IntegerMatrix *));
+ if ( (win == NULL) || (ser == NULL) || (mat == NULL) ) {
ERROR("Failed to expand series buffers\n");
return 1;
}
@@ -237,7 +330,8 @@ int main(int argc, char *argv[])
close_stream(st);
free(win);
- free(series);
+ free(ser);
+ free(mat);
return 0;
}