From d9153e1c7d5cd0fb46211aad627020341c13f07d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 17 Oct 2014 17:20:37 +0200 Subject: Find rotation series even when images arrive out of order --- src/whirligig.c | 238 ++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 143 insertions(+), 95 deletions(-) (limited to 'src/whirligig.c') diff --git a/src/whirligig.c b/src/whirligig.c index 4b968027..4a76d1b9 100644 --- a/src/whirligig.c +++ b/src/whirligig.c @@ -246,45 +246,50 @@ static int gatinator(UnitCell *a, UnitCell *b, IntegerMatrix **pmb) } -/* Try all combinations of crystals from the two patterns, in the hope of - * finding the start of a rotation series */ -static int try_all(struct image *a, struct image *b, int *c1, int *c2, - IntegerMatrix **m2) +static int try_all(struct image *win, signed int *ser, IntegerMatrix **m, + int ws, signed int p1, signed int p2) { int i, j; - for ( i=0; in_crystals; i++ ) { - for ( j=0; jn_crystals; j++ ) { - if ( gatinator(crystal_get_cell(a->crystals[i]), - crystal_get_cell(b->crystals[j]), m2) ) - { - *c1 = i; - *c2 = j; - return 1; - } + for ( i=0; icrystals[c1]), tfn); + if ( ser[p2] == -1 ) { + return try_all(win, ser, m, ws, p1, p2); + } + + tfn = tfn_from_intmat(m[p2]); + ref = cell_transform(crystal_get_cell(win[p2].crystals[ser[p2]]), tfn); tfn_free(tfn); - for ( j=0; jn_crystals; j++ ) { - if ( gatinator(ref, crystal_get_cell(b->crystals[j]), m2) ) { - *c2 = j; + for ( j=0; j 0) && (win[pos-1].serial != 0) ) { + try_join(win, ser, m, ws, pos, pos-1); + } + + /* Try to connect to the right */ + if ( (pos < ws-1) && (win[pos+1].serial != 0) ) { + try_join(win, ser, m, ws, pos, pos+1); + } +} + + static void dump(struct image *win, signed int *ser, IntegerMatrix **mat, int window_len, int pos) { @@ -312,6 +333,78 @@ static void dump(struct image *win, signed int *ser, IntegerMatrix **mat, } +static int series_fills_window(signed int *ser, int ws) +{ + int i; + + for ( i=0; iserial > win[ws-1].serial ) { + + int i; + + if ( series_fills_window(ser, ws) ) { + + ws++; + win = realloc(win, ws*sizeof(struct image)); + ser = realloc(ser, ws*sizeof(signed int)); + mat = realloc(mat, ws*sizeof(IntegerMatrix *)); + if ( (win == NULL) || (ser == NULL) || (mat == NULL) ) { + ERROR("Failed to expand series buffers\n"); + exit(1); + } + *pws = ws; + *pwin = win; + *pser = ser; + *pmat = mat; + + pos = ws-1; + + } else { + + sf = cur->serial - win[ws-1].serial; + + memmove(win, win+sf, (ws-sf)*sizeof(struct image)); + memmove(ser, ser+sf, (ws-sf)*sizeof(signed int)); + memmove(mat, mat+sf, (ws-sf)*sizeof(IntegerMatrix *)); + + for ( i=0; iserial)-1; + + } + + win[pos] = *cur; + + return pos; +} + + static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); @@ -329,11 +422,11 @@ int main(int argc, char *argv[]) int c; Stream *st; int polarisation = 1; - int pos = 0; struct image *win; signed int *ser; IntegerMatrix **mat; - int window_len = 64; + int default_window_size = 16; + int ws, i; /* Long options */ const struct option longopts[] = { @@ -385,102 +478,57 @@ int main(int argc, char *argv[]) return 1; } - win = calloc(window_len, sizeof(struct image)); - ser = calloc(window_len, sizeof(int)); - mat = calloc(window_len, sizeof(IntegerMatrix *)); + /* Allocate initial window */ + ws = default_window_size; + win = calloc(ws, sizeof(struct image)); + ser = calloc(ws, sizeof(signed int)); + mat = calloc(ws, sizeof(IntegerMatrix *)); if ( (win == NULL) || (ser == NULL) || (mat == NULL) ) { ERROR("Failed to allocate series buffers\n"); return 1; } - pos = 0; - do { + for ( i=0; idiv = NAN; - cur->bw = NAN; - cur->det = NULL; - if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS + cur.div = NAN; + cur.bw = NAN; + cur.det = NULL; + if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS | STREAM_READ_UNITCELL) != 0 ) { break; } - if ( isnan(cur->div) || isnan(cur->bw) ) { + if ( isnan(cur.div) || isnan(cur.bw) ) { ERROR("Chunk doesn't contain beam parameters.\n"); return 1; } - /* Need at least two images to compare */ - if ( pos == 0 ) { - ser[0] = -1; - pos++; - continue; - } - - if ( ser[pos-1] == -1 ) { - if ( try_all(&win[pos-1], cur, &c1, &c2, &mat[pos]) ) { - ser[pos-1] = c1; - ser[pos] = c2; - } else { - ser[pos] = -1; - mat[pos] = NULL; - } - } else { - if ( try_cont(&win[pos-1], cur, ser[pos-1], &c2, - mat[pos-1], &mat[pos]) ) - { - ser[pos] = c2; - } else { - ser[pos] = -1; - mat[pos] = NULL; - } - } - - if ( ser[pos-1] != -1 ) { - printf("*"); - } else { - printf("-"); - } - fflush(stdout); - - if ( ser[0] == -1 ) { - dump(win, ser, mat, window_len, 1); - pos--; + if ( cur.serial < 1 ) { + ERROR("Serial numbers must be greater than zero.\n"); + return 1; } - if ( (pos > 0) && (ser[pos] == -1) && (ser[pos-1] == -1) ) { - /* Series ready to process */ - process_series(win, ser, mat, pos-1); - dump(win, ser, mat, window_len, pos); - pos = 0; - } - - pos++; - if ( pos == window_len ) { - window_len *= 2; - win = realloc(win, window_len*sizeof(struct image)); - 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; - } - } + pos = add_to_window(&cur, &win, &ser, &mat, &ws); + try_connect(win, ser, mat, ws, pos); } while ( 1 ); printf("\n"); close_stream(st); - /* Final series to process? */ - if ( pos > 0 ) { - process_series(win, ser, mat, pos); - dump(win, ser, mat, window_len, 1); + for ( i=0; i