diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/Makefile.am | 5 | ||||
-rw-r--r-- | src/basis.c | 12 | ||||
-rw-r--r-- | src/cache.c | 17 | ||||
-rw-r--r-- | src/control.c | 1 | ||||
-rw-r--r-- | src/control.h | 22 | ||||
-rw-r--r-- | src/dirax.c | 437 | ||||
-rw-r--r-- | src/dirax.h | 28 | ||||
-rw-r--r-- | src/displaywindow.c | 119 | ||||
-rw-r--r-- | src/ipr.c | 2 | ||||
-rw-r--r-- | src/main.c | 30 | ||||
-rw-r--r-- | src/reflections.c | 20 | ||||
-rw-r--r-- | src/reflections.h | 21 | ||||
-rw-r--r-- | src/utils.c | 28 | ||||
-rw-r--r-- | src/utils.h | 3 |
14 files changed, 641 insertions, 104 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index bdfe197..0be30fc 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,7 +1,8 @@ bin_PROGRAMS = dtr dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c cache.c \ - itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c -dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas + itrans-threshold.c itrans-zaefferer.c itrans-lsq.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \ + dirax.c +dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" diff --git a/src/basis.c b/src/basis.c index 19004ee..8a282e3 100644 --- a/src/basis.c +++ b/src/basis.c @@ -87,7 +87,7 @@ static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) { Reflection *check; - check = reflection_find_nearest(ctx->reflectionlist, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j); + check = reflectionlist_find_nearest(ctx->reflectionlist, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j); if ( check && (distance3d(check->x, check->y, check->z, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j) < tol) ) { lfom++; } @@ -135,14 +135,14 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { z = z_temp; /* Find the point in the middle of the "wedge" */ - scale = reflection_largest_g(ctx->reflectionlist)/6; + scale = reflectionlist_largest_g(ctx->reflectionlist)/6; x *= scale; y *= scale; z *= scale; reflection_add(ctx->reflectionlist, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2); /* Find an "origin" reflection */ - centre = reflection_find_nearest(ctx->reflectionlist, x, y, z); + centre = reflectionlist_find_nearest(ctx->reflectionlist, x, y, z); if ( !centre ) return NULL; centre->found = 1; reflection_add(ctx->reflectionlist, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED); @@ -161,7 +161,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { accept = 1; /* Find a "candidate vector" reflection */ - vector = reflection_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */ + vector = reflectionlist_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); if ( !vector ) { printf("BS: Only found %i seeds\n", i); return seeds; @@ -174,7 +174,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { vz = vector->z - centre->z; /* Proximity test: don't duplicate seeds */ - check = reflection_find_nearest(seeds, vx, vy, vz); + check = reflectionlist_find_nearest(seeds, vx, vy, vz); if ( check ) { if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) { /* Too close to another seed */ @@ -209,7 +209,7 @@ static Basis *basis_assemble_seeds(ReflectionList *seeds) { int i, j, b; ReflectionList *seeds_sorted; - seeds_sorted = reflection_sort_lfom(seeds); + seeds_sorted = reflectionlist_sort_lfom(seeds); basis = malloc(10*sizeof(Basis)); diff --git a/src/cache.c b/src/cache.c index 078fa72..ad47ca2 100644 --- a/src/cache.c +++ b/src/cache.c @@ -41,10 +41,10 @@ ReflectionList *cache_load(const char *filename) { reflectionlist = reflectionlist_new(); f = fopen(filename, "rb"); if ( !f ) { - fprintf(stderr, "Couldn't read reflections from cache\n"); + fprintf(stderr, "Couldn't open cache file\n"); } if ( fread(&ch, sizeof(CacheHeader), 1, f) == 0 ) { - fprintf(stderr, "Couldn't read reflections from cache\n"); + fprintf(stderr, "Couldn't read cache header\n"); fclose(f); return NULL; } @@ -58,7 +58,7 @@ ReflectionList *cache_load(const char *filename) { fprintf(stderr, "Couldn't read reflections from cache\n"); fclose(f); free(cr); - reflection_clear(reflectionlist); + reflectionlist_clear(reflectionlist); return NULL; } @@ -129,9 +129,13 @@ unsigned int cache_is_cachefile(const char *filename) { size_t nread; fh = fopen(filename, "rb"); + if ( !fh ) { + printf("Couldn't open file '%s'\n", filename); + return 0; + } nread = fread(&ch, sizeof(CacheHeader), 1, fh); fclose(fh); - + if ( nread != 1 ) { return 0; } @@ -140,11 +144,6 @@ unsigned int cache_is_cachefile(const char *filename) { return 1; } - /* Backwards compatability */ - if ( strncmp(ch.top, "DTR-CACHE", 9) == 0 ) { - return 1; - } - return 0; } diff --git a/src/control.c b/src/control.c index 270be6f..032700d 100644 --- a/src/control.c +++ b/src/control.c @@ -53,6 +53,7 @@ ControlContext *control_ctx_new() { ctx->x_centre = 0; ctx->y_centre = 0; ctx->have_centres = 0; + ctx->cell = NULL; return ctx; diff --git a/src/control.h b/src/control.h index 5b0d187..bb42395 100644 --- a/src/control.h +++ b/src/control.h @@ -26,7 +26,8 @@ typedef enum ift_enum { INPUT_NONE, INPUT_QDRP, INPUT_MRC, - INPUT_CACHE + INPUT_CACHE, + INPUT_DRX } InputFileType; typedef enum { @@ -44,8 +45,9 @@ typedef enum { } PeakSearchMode; typedef enum { - RECONSTRUCTION_MAPPING, - RECONSTRUCTION_PREDICTION + RECONSTRUCTION_MAPPING = 0, + RECONSTRUCTION_PREDICTION = 1<<0, + RECONSTRUCTION_DIRAX = 1<<1 } ReconstructionMode; typedef struct imagerecord_struct { @@ -85,6 +87,7 @@ typedef struct cctx_struct { unsigned int prealign; unsigned int savecache; unsigned int have_centres; + unsigned int use_dirax; /* Input filename */ char *filename; @@ -113,18 +116,31 @@ typedef struct cctx_struct { /* Output */ struct reflectionlist_struct *reflectionlist; struct dw_struct *dw; + struct basis_struct *cell; /* GTK bits */ GtkWidget *combo_algorithm; GtkWidget *combo_peaksearch; GtkWidget *checkbox_prealign; GtkWidget *checkbox_savecache; + GtkWidget *checkbox_dirax; /* IPR stuff */ int ipr_cur_image; struct imagedisplay_struct *ipr_id; struct reflectionlist_struct *ipr_lat; struct basis_struct *ipr_basis; + + /* DirAx low-level stuff */ + GIOChannel *dirax; + int dirax_pty; + char *dirax_rbuffer; + int dirax_rbufpos; + int dirax_rbuflen; + + /* DirAx high-level stuff */ + int dirax_step; + int dirax_read_cell; } ControlContext; diff --git a/src/dirax.c b/src/dirax.c new file mode 100644 index 0000000..4f4d135 --- /dev/null +++ b/src/dirax.c @@ -0,0 +1,437 @@ +/* + * dirax.c + * + * Invoke the DirAx auto-indexing program + * also: handle DirAx input files + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <glib.h> +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <string.h> +#include <pty.h> +#include <unistd.h> +#include <sys/wait.h> +#include <fcntl.h> +#include <assert.h> +#include <sys/ioctl.h> +#include <termio.h> + +#include "control.h" +#include "reflections.h" +#include "utils.h" +#include "basis.h" +#include "displaywindow.h" + +typedef enum { + DIRAX_INPUT_NONE, + DIRAX_INPUT_LINE, + DIRAX_INPUT_PROMPT +} DirAxInputType; + +static void dirax_parseline(const char *line, ControlContext *ctx) { + + int i, rf; + char *copy; + + copy = strdup(line); + for ( i=0; i<strlen(copy); i++ ) { + if ( copy[i] == '\r' ) copy[i]='r'; + if ( copy[i] == '\n' ) copy[i]='n'; + } + printf("DX: DirAx: %s\n", copy); + free(copy); + + /* Is this the first line of a unit cell specification? */ + rf = 0; i = 0; + while ( (i<strlen(line)) && ((line[i] == 'R') || (line[i] == 'D') || (line[i] == ' ')) ) { + if ( line[i] == 'R' ) rf = 1; + if ( (line[i] == 'D') && rf ) { + ctx->dirax_read_cell = 1; + if ( ctx->cell ) { + free(ctx->cell); + } + ctx->cell = malloc(sizeof(Basis)); + ctx->cell->a.x = 0.0; ctx->cell->a.y = 0.0; ctx->cell->a.z = 0.0; + ctx->cell->b.x = 0.0; ctx->cell->b.y = 0.0; ctx->cell->b.z = 0.0; + ctx->cell->c.x = 0.0; ctx->cell->c.y = 0.0; ctx->cell->c.z = 0.0; + return; + } + i++; + } + + /* Parse unit cell vectors as appropriate */ + if ( ctx->dirax_read_cell == 1 ) { + /* First row of unit cell values */ + float x, y, z; + sscanf(line, "%f %f %f", &x, &y, &z); + ctx->cell->a.x = x*1e10; ctx->cell->a.y = y*1e10; ctx->cell->a.z = z*1e10; + ctx->dirax_read_cell++; + return; + } else if ( ctx->dirax_read_cell == 2 ) { + /* First row of unit cell values */ + float x, y, z; + sscanf(line, "%f %f %f", &x, &y, &z); + ctx->cell->b.x = x*1e10; ctx->cell->b.y = y*1e10; ctx->cell->b.z = z*1e10; + ctx->dirax_read_cell++; + return; + } else if ( ctx->dirax_read_cell == 3 ) { + /* First row of unit cell values */ + float x, y, z; + sscanf(line, "%f %f %f", &x, &y, &z); + ctx->cell->c.x = x*1e10; ctx->cell->c.y = y*1e10; ctx->cell->c.z = z*1e10; + printf("DX: Read a reciprocal unit cell\n"); + displaywindow_update(ctx->dw); + ctx->dirax_read_cell = 0; + return; + } + + ctx->dirax_read_cell = 0; + +} + +static void dirax_sendline(const char *line, ControlContext *ctx) { + + char *copy; + int i; + + write(ctx->dirax_pty, line, strlen(line)); + + copy = strdup(line); + for ( i=0; i<strlen(copy); i++ ) { + if ( copy[i] == '\r' ) copy[i]='\0'; + if ( copy[i] == '\n' ) copy[i]='\0'; + } + printf("sent '%s'", copy); /* No newline here */ + free(copy); + +} + +static void dirax_send_next(ControlContext *ctx) { + + switch ( ctx->dirax_step ) { + + case 0 : { + dirax_sendline("read dtr.drx\n", ctx); + ctx->dirax_step++; + break; + } + + case 1 : { + dirax_sendline("go\n", ctx); + ctx->dirax_step++; + break; + } + + case 2 : { + dirax_sendline("cell\n", ctx); + ctx->dirax_step++; + break; + } + + default: { + printf("idle"); + } + + } + + fsync(ctx->dirax_pty); + +} + +static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, ControlContext *ctx) { + + int rval; + + rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos, ctx->dirax_rbuflen-ctx->dirax_rbufpos); + + if ( (rval == -1) || (rval == 0) ) { + + printf("DX: Lost connection to DirAx\n"); + g_io_channel_shutdown(ctx->dirax, FALSE, NULL); + ctx->dirax = NULL; + return FALSE; + + } else { + + int no_string = 0; + + ctx->dirax_rbufpos += rval; + assert(ctx->dirax_rbufpos <= ctx->dirax_rbuflen); + + while ( (!no_string) && (ctx->dirax_rbufpos > 0) ) { + + int i; + int block_ready = 0; + DirAxInputType type = DIRAX_INPUT_NONE; + + /* See if there's a full line in the buffer yet */ + for ( i=0; i<ctx->dirax_rbufpos; i++ ) { /* Means the last value looked at is roffset-1 */ + + if ( (ctx->dirax_rbuffer[i] == '\n') ) { + block_ready = 1; + type = DIRAX_INPUT_LINE; + break; + } + + if ( i <= ctx->dirax_rbufpos-7 ) { + if ( strncmp(ctx->dirax_rbuffer+i, "Dirax> ", 7) == 0 ) { + block_ready = 1; + type = DIRAX_INPUT_PROMPT; + break; + } + } + + } + + if ( block_ready == 1 ) { + + switch ( type ) { + + case DIRAX_INPUT_LINE : { + + char *block_buffer = NULL; + unsigned int new_rbuflen; + unsigned int endbit_length; + + block_buffer = malloc(i+1); + memcpy(block_buffer, ctx->dirax_rbuffer, i); + block_buffer[i] = '\0'; + dirax_parseline(block_buffer, ctx); + free(block_buffer); + + /* Now the block's been parsed, it should be forgotten about */ + endbit_length = i+1; + memmove(ctx->dirax_rbuffer, ctx->dirax_rbuffer + endbit_length, + ctx->dirax_rbuflen - endbit_length); + /* Subtract the number of bytes removed */ + ctx->dirax_rbufpos = ctx->dirax_rbufpos - endbit_length; + new_rbuflen = ctx->dirax_rbuflen - endbit_length; + if ( new_rbuflen == 0 ) { + new_rbuflen = 256; + } + ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, new_rbuflen); + ctx->dirax_rbuflen = new_rbuflen; + + break; + + } + + case DIRAX_INPUT_PROMPT : { + + memmove(ctx->dirax_rbuffer+i, ctx->dirax_rbuffer+i+7, 7); + ctx->dirax_rbufpos -= 7; + printf("DX: Got prompt: "); + dirax_send_next(ctx); + printf("\n"); + + break; + + } + + default : { + printf("DX: Unrecognised input mode (this never happens!)\n"); + abort(); + } + + } + + } else { + + if ( ctx->dirax_rbufpos == ctx->dirax_rbuflen ) { + + /* More buffer space is needed */ + ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, ctx->dirax_rbuflen + 256); + ctx->dirax_rbuflen = ctx->dirax_rbuflen + 256; + /* The new space gets used at the next read, shortly... */ + + } + no_string = 1; + + } + + } + + } + + return TRUE; + +} + +void dirax_invoke(ControlContext *ctx) { + + FILE *fh; + Reflection *ref; + int pty; + unsigned int opts; + pid_t dirax_pid; + + printf("DX: Starting DirAx...\n"); + + fh = fopen("dtr.drx", "w"); + if ( !fh ) { + printf("DX: Couldn't open temporary file dtr.drx\n"); + return; + } + fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. DirAx can't handle the truth. */ + ref = ctx->reflectionlist->reflections; + while ( ref ) { + fprintf(fh, "%f %f %f %f\n", ref->x/1e10, ref->y/1e10, ref->z/1e10, ref->intensity); + ref = ref->next; + } + fclose(fh); + + struct termios tty; + if ( ioctl(0, TCGETS, &tty) ) { + printf("DX: Failed to read pty attributes\n"); + return; + } + tty.c_lflag &= ~ECHO; + dirax_pid = forkpty(&pty, NULL, &tty, NULL); + if ( dirax_pid == -1 ) { + printf("DX: Failed to fork.\n"); + return; + } + if ( dirax_pid == 0 ) { + + /* Child process: invoke DirAx */ + execlp("dirax", "", (char *)NULL); + printf("(from the Other Side) Failed to invoke DirAx.\n"); + _exit(0); + + } + + ctx->dirax_pty = pty; + ctx->dirax_rbuffer = malloc(256); + ctx->dirax_rbuflen = 256; + ctx->dirax_rbufpos = 0; + + /* Set non-blocking */ + opts = fcntl(pty, F_GETFL); + fcntl(pty, F_SETFL, opts | O_NONBLOCK); + + ctx->dirax_step = 0; + ctx->dirax_read_cell = 0; + ctx->dirax = g_io_channel_unix_new(pty); + g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable, ctx); + + return; + +} + +/* Despite being part of the same module, this has very little to do with invoking DirAx */ +ReflectionList *dirax_load(const char *filename) { + + FILE *fh; + char line[256]; + ReflectionList *list; + int lambda_set = 0; + + fh = fopen(filename, "r"); + if ( !fh ) { + printf("Couldn't open file '%s'\n", filename); + return 0; + } + + list = reflectionlist_new(); + + while ( !feof(fh) ) { + + size_t ptr; + float lambda, theta, chib, phib; + + fgets(line, 255, fh); + ptr = skipspace(line); + if ( line[ptr] == '!' ) continue; + if ( line[ptr] == '\n' ) continue; + if ( line[ptr] == '\r' ) continue; + if ( sscanf(line+ptr, "%f %f %f\n", &theta, &phib, &chib) == 3 ) { + + double s, x, y, z; + float blah, intensity; + + /* Try to find an intensity value. Use dummy value if it fails */ + if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah, &blah, &intensity) != 4 ) { + intensity = 1.0; + } + + if ( !lambda_set ) { + printf("DX: Wavelength not specified\n"); + continue; + } + + chib = deg2rad(chib); + phib = deg2rad(phib); + theta = deg2rad(theta); + s = 2*(sin(theta)/lambda); + x = -s*cos(chib)*sin(phib); + y = +s*cos(chib)*cos(phib); + z = +s*sin(chib); + reflection_add(list, x, y, z, 1.0, REFLECTION_NORMAL); + + continue; + + } + if ( sscanf(line+ptr, "%f\n", &lambda) == 1 ) { + if ( lambda_set ) { + printf("DX: Warning: Found something which looks like a second wavelength\n"); + } + lambda /= 1e10; /* Convert from A to m */ + lambda_set = 1; + } + + } + + fclose(fh); + + return list; + +} + +int dirax_is_drxfile(const char *filename) { + + FILE *fh; + float lambda; + char line[256]; + + fh = fopen(filename, "r"); + if ( !fh ) { + printf("Couldn't open file '%s'\n", filename); + return 0; + } + + while ( !feof(fh) ) { + + size_t ptr; + + fgets(line, 255, fh); + ptr = skipspace(line); + if ( line[ptr] == '!' ) continue; + if ( line[ptr] == '\n' ) continue; + if ( line[ptr] == '\r' ) continue; + fscanf(fh, "%f\n", &lambda); + fclose(fh); + if ( lambda > 0.5 ) { + return 1; + } else { + return 0; + } + + } + + fclose(fh); + + return 0; + +} + diff --git a/src/dirax.h b/src/dirax.h new file mode 100644 index 0000000..b06c6c9 --- /dev/null +++ b/src/dirax.h @@ -0,0 +1,28 @@ +/* + * dirax.h + * + * Invoke the DirAx auto-indexing program + * also: handle DirAx input files + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifndef DIRAX_H +#define DIRAX_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "control.h" +#include "reflections.h" + +extern void dirax_invoke(ControlContext *ctx); +extern ReflectionList *dirax_load(const char *filename); +extern int dirax_is_drxfile(const char *filename); + +#endif /* DIRAX_H */ + diff --git a/src/displaywindow.c b/src/displaywindow.c index b38007d..866f63f 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -31,13 +31,14 @@ #include "ipr.h" #include "displaywindow.h" #include "basis.h" +#include "dirax.h" enum { DW_ORTHO, DW_PERSPECTIVE }; -static void displaywindow_about() { +static void displaywindow_about(GtkWidget *widget, DisplayWindow *dw) { GtkWidget *window; @@ -70,10 +71,14 @@ static void displaywindow_about() { } -static void displaywindow_close() { +static void displaywindow_close(GtkWidget *widget, DisplayWindow *dw) { gtk_exit(0); } +static void displaywindow_dirax(GtkWidget *widget, DisplayWindow *dw) { + dirax_invoke(dw->ctx); +} + static void displaywindow_gl_set_ortho(DisplayWindow *dw, GLfloat w, GLfloat h) { glMatrixMode(GL_PROJECTION); @@ -169,6 +174,7 @@ static void displaywindow_addmenubar(DisplayWindow *dw) { { "ViewAction", NULL, "_View", NULL, NULL, NULL }, { "ToolsAction", NULL, "_Tools", NULL, NULL, NULL }, + { "DirAxAction", GTK_STOCK_EXECUTE, "Invoke _DirAx", "<Ctrl>D", NULL, G_CALLBACK(displaywindow_dirax) }, { "HelpAction", NULL, "_Help", NULL, NULL, NULL }, { "AboutAction", GTK_STOCK_ABOUT, "_About DTR...", NULL, NULL, G_CALLBACK(displaywindow_about) }, @@ -343,7 +349,6 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { GLfloat green[] = { 0.0, 1.0, 0.0, 1.0 }; GLfloat yellow[] = { 1.0, 1.0, 0.0, 1.0 }; GLfloat glass[] = { 0.2, 0.0, 0.8, 000.1 }; - GLfloat glass_spec[] = { 0.8, 0.8, 0.8, 1.0 }; GLfloat black[] = { 0.0, 0.0, 0.0, 1.0 }; Reflection *reflection; GLUquadricObj *quadric; @@ -360,15 +365,15 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { } reflection = ctx->reflectionlist->reflections; i = 0; - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_NORMAL ) i++; reflection = reflection->next; - } while ( reflection != NULL ); + }; dw->gl_ref_num_vertices = i; i = 0; reflection = ctx->reflectionlist->reflections; vertices = malloc(3*dw->gl_ref_num_vertices*sizeof(GLfloat)); - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_NORMAL ) { vertices[3*i + 0] = reflection->x/1e9; vertices[3*i + 1] = reflection->y/1e9; @@ -376,7 +381,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { i++; } reflection = reflection->next; - } while ( reflection != NULL ); + }; if ( dw->gl_use_buffers ) { glBindBufferARB(GL_ARRAY_BUFFER, dw->gl_ref_vertex_buffer); glBufferDataARB(GL_ARRAY_BUFFER, 3*dw->gl_ref_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW); @@ -392,21 +397,21 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { } reflection = ctx->reflectionlist->reflections; i = 0; - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_MARKER ) i++; reflection = reflection->next; - } while ( reflection != NULL ); + }; dw->gl_marker_num_vertices = i*VERTICES_IN_A_BLOB; i = 0; reflection = ctx->reflectionlist->reflections; vertices = malloc(3*dw->gl_marker_num_vertices*sizeof(GLfloat)); normals = malloc(3*dw->gl_marker_num_vertices*sizeof(GLfloat)); - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_MARKER ) { DRAW_BLOB } reflection = reflection->next; - } while ( reflection != NULL ); + }; if ( dw->gl_use_buffers ) { glBindBufferARB(GL_ARRAY_BUFFER, dw->gl_marker_vertex_buffer); glBufferDataARB(GL_ARRAY_BUFFER, 3*dw->gl_marker_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW); @@ -426,21 +431,21 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { } reflection = ctx->reflectionlist->reflections; i = 0; - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_GENERATED ) i++; reflection = reflection->next; - } while ( reflection != NULL ); + }; dw->gl_gen_num_vertices = i*VERTICES_IN_A_BLOB; i = 0; reflection = ctx->reflectionlist->reflections; vertices = malloc(3*dw->gl_gen_num_vertices*sizeof(GLfloat)); normals = malloc(3*dw->gl_gen_num_vertices*sizeof(GLfloat)); - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_GENERATED ) { DRAW_BLOB } reflection = reflection->next; - } while ( reflection != NULL ); + }; if ( dw->gl_use_buffers ) { glBindBufferARB(GL_ARRAY_BUFFER, dw->gl_gen_vertex_buffer); glBufferDataARB(GL_ARRAY_BUFFER, 3*dw->gl_gen_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW); @@ -519,7 +524,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { /* Plot the other reflections */ reflection = ctx->reflectionlist->reflections; quadric = gluNewQuadric(); - do { + while ( reflection != NULL ) { if ( reflection->type == REFLECTION_CENTRAL ) { @@ -583,61 +588,61 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { reflection = reflection->next; - } while ( reflection != NULL ); + }; - /* If this is an iterative prediction-refinement reconstruction, draw the unit cell */ - if ( ctx->rmode == RECONSTRUCTION_PREDICTION ) { + /* Draw the reciprocal unit cell if one is available */ + if ( ctx->cell ) { glBegin(GL_LINE_STRIP); glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black); glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0); glVertex3f(0.0, 0.0, 0.0); - glVertex3f(ctx->ipr_basis->a.x/1e9, ctx->ipr_basis->a.y/1e9, ctx->ipr_basis->a.z/1e9); - glVertex3f(ctx->ipr_basis->a.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->a.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->a.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->b.x/1e9, ctx->ipr_basis->b.y/1e9, ctx->ipr_basis->b.z/1e9); + glVertex3f(ctx->cell->a.x/1e9, ctx->cell->a.y/1e9, ctx->cell->a.z/1e9); + glVertex3f(ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->b.x/1e9, ctx->cell->b.y/1e9, ctx->cell->b.z/1e9); glVertex3f(0.0, 0.0, 0.0); - glVertex3f(ctx->ipr_basis->c.x/1e9, ctx->ipr_basis->c.y/1e9, ctx->ipr_basis->c.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->a.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->a.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->a.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->a.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->a.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->a.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9, ctx->ipr_basis->c.y/1e9, ctx->ipr_basis->c.z/1e9); + glVertex3f(ctx->cell->c.x/1e9, ctx->cell->c.y/1e9, ctx->cell->c.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9, ctx->cell->c.y/1e9, ctx->cell->c.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->b.x/1e9, ctx->ipr_basis->b.y/1e9, ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->a.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->a.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->a.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->a.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->a.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->a.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->a.x/1e9 + ctx->ipr_basis->b.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->a.y/1e9 + ctx->ipr_basis->b.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->a.z/1e9 + ctx->ipr_basis->b.z/1e9); - glVertex3f(ctx->ipr_basis->c.x/1e9 + ctx->ipr_basis->a.x/1e9, - ctx->ipr_basis->c.y/1e9 + ctx->ipr_basis->a.y/1e9, - ctx->ipr_basis->c.z/1e9 + ctx->ipr_basis->a.z/1e9); - glVertex3f(ctx->ipr_basis->a.x/1e9, ctx->ipr_basis->a.y/1e9, ctx->ipr_basis->a.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->b.x/1e9, ctx->cell->b.y/1e9, ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9); + glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9, + ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9, + ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9); + glVertex3f(ctx->cell->a.x/1e9, ctx->cell->a.y/1e9, ctx->cell->a.z/1e9); glEnd(); } /* Zero plane */ glBegin(GL_QUADS); glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, glass); - glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, glass_spec); - glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 50.0); + glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black); + glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0); glVertex3f(50, 50, 0.0); glVertex3f(50, -50, 0.0); glVertex3f(-50, -50, 0.0); @@ -95,7 +95,7 @@ static gint ipr_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext if ( ctx->ipr_cur_image == ctx->n_images ) ctx->ipr_cur_image = 0; imagedisplay_clear_circles(ctx->ipr_id); - reflection_clear_markers(ctx->reflectionlist); + reflectionlist_clear_markers(ctx->reflectionlist); refl = reproject_get_reflections(ctx->images[ctx->ipr_cur_image], &n, ctx->ipr_lat, ctx); for ( j=0; j<n; j++ ) { @@ -33,12 +33,13 @@ #include "mapping.h" #include "prealign.h" #include "control.h" +#include "dirax.h" void main_do_reconstruction(ControlContext *ctx) { int val = 0; - if ( ctx->inputfiletype != INPUT_CACHE ) { + if ( (ctx->inputfiletype != INPUT_CACHE) && (ctx->inputfiletype != INPUT_DRX) ) { prealign_sum_stack(ctx); mapping_create(ctx); } @@ -64,7 +65,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response, if ( response == GTK_RESPONSE_OK ) { - int val = -1; + int val = 0; switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_algorithm)) ) { case 0 : ctx->rmode = RECONSTRUCTION_MAPPING; break; @@ -103,7 +104,10 @@ static gint main_method_window_response(GtkWidget *method_window, gint response, val = mrc_read(ctx); } else if ( ctx->inputfiletype == INPUT_CACHE ) { ctx->reflectionlist = cache_load(ctx->filename); - val=0; + if ( !ctx->reflectionlist ) val = 1; + } else if ( ctx->inputfiletype == INPUT_DRX ) { + ctx->reflectionlist = dirax_load(ctx->filename); + if ( !ctx->reflectionlist ) val = 1; } if ( val ) { @@ -174,7 +178,12 @@ void main_method_dialog_open(ControlContext *ctx) { gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_savecache, 1, 3, 4, 5); if ( ctx->inputfiletype == INPUT_CACHE ) { - gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get from cache file"); + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "3D coordinates from cache file"); + gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch), FALSE); + gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 5); + } + if ( ctx->inputfiletype == INPUT_DRX ) { + gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "3D coordinates from DirAx file"); gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch), FALSE); gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 5); } @@ -194,6 +203,7 @@ int main(int argc, char *argv[]) { ControlContext *ctx; InputFileType type; struct stat stat_buffer; + FILE *fh; gtk_init(&argc, &argv); @@ -212,10 +222,17 @@ int main(int argc, char *argv[]) { return 1; } - filename = basename(argv[1]); + filename = argv[1]; ctx = control_ctx_new(); type = INPUT_NONE; + fh = fopen(filename, "r"); + if ( !fh ) { + printf("Couldn't open file '%s'\n", filename); + return 1; + } + fclose(fh); + if ( qdrp_is_qdrprc(filename) ) { printf("QDRP input file detected.\n"); ctx->inputfiletype = INPUT_QDRP; @@ -225,6 +242,9 @@ int main(int argc, char *argv[]) { } else if ( cache_is_cachefile(filename) ) { printf("Cached reflection file detected.\n"); ctx->inputfiletype = INPUT_CACHE; + } else if ( dirax_is_drxfile(filename) ) { + printf("Dirax input file detected.\n"); + ctx->inputfiletype = INPUT_DRX; } else { fprintf(stderr, "Unrecognised input file type\n"); return 1; diff --git a/src/reflections.c b/src/reflections.c index 56871f0..dc3ce2f 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -39,7 +39,7 @@ ReflectionList *reflectionlist_new() { } -void reflection_clear_markers(ReflectionList *reflectionlist) { +void reflectionlist_clear_markers(ReflectionList *reflectionlist) { Reflection *reflection = reflectionlist->reflections; Reflection *prev = NULL; @@ -70,7 +70,7 @@ void reflection_clear_markers(ReflectionList *reflectionlist) { reflectionlist->last_reflection = prev; } -void reflection_clear(ReflectionList *reflectionlist) { +void reflectionlist_clear(ReflectionList *reflectionlist) { Reflection *reflection = reflectionlist->reflections; do { @@ -83,8 +83,8 @@ void reflection_clear(ReflectionList *reflectionlist) { } -void reflection_free(ReflectionList *reflectionlist) { - reflection_clear(reflectionlist); +void reflectionlist_free(ReflectionList *reflectionlist) { + reflectionlist_clear(reflectionlist); free(reflectionlist); } @@ -101,7 +101,7 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d reflectionlist->list_capped = 1; } - nearest = reflection_find_nearest_type(reflectionlist, x, y, z, type); + nearest = reflectionlist_find_nearest_type(reflectionlist, x, y, z, type); if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL; new_reflection = malloc(sizeof(Reflection)); @@ -226,7 +226,7 @@ void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection * reflectionlist->n_reflections++; } -double reflection_largest_g(ReflectionList *reflectionlist) { +double reflectionlist_largest_g(ReflectionList *reflectionlist) { double max = 0.0; Reflection *reflection; @@ -245,7 +245,7 @@ double reflection_largest_g(ReflectionList *reflectionlist) { } -Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, double y, double z) { +Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x, double y, double z) { double max = +INFINITY; Reflection *reflection; @@ -268,7 +268,7 @@ Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, do } -Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) { +Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) { double max = +INFINITY; Reflection *reflection; @@ -291,7 +291,7 @@ Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlis } -Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type) { +Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type) { double max = +INFINITY; Reflection *reflection; @@ -316,7 +316,7 @@ Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double } /* This destroys the lfom values in the input list */ -ReflectionList *reflection_sort_lfom(ReflectionList *in) { +ReflectionList *reflectionlist_sort_lfom(ReflectionList *in) { ReflectionList *out; Reflection *best; diff --git a/src/reflections.h b/src/reflections.h index c753430..11fbdcf 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -59,22 +59,21 @@ typedef struct reflectionlist_struct { } ReflectionList; extern ReflectionList *reflectionlist_new(void); -extern void reflection_clear(ReflectionList *reflectionlist); -extern void reflection_clear_markers(ReflectionList *reflectionlist); -extern void reflection_free(ReflectionList *reflectionlist); - -extern Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, double z, double intensity, ReflectionType type); +extern void reflectionlist_clear(ReflectionList *reflectionlist); +extern void reflectionlist_clear_markers(ReflectionList *reflectionlist); +extern void reflectionlist_free(ReflectionList *reflectionlist); #include "control.h" - +extern Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, double z, double intensity, ReflectionType type); extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord *imagerecord, double intensity); extern void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r); -extern double reflection_largest_g(ReflectionList *reflectionlist); extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta); -extern Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, double y, double z); -extern Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance); -extern Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type); -extern ReflectionList *reflection_sort_lfom(ReflectionList *in); + +extern double reflectionlist_largest_g(ReflectionList *reflectionlist); +extern Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x, double y, double z); +extern Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance); +extern Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type); +extern ReflectionList *reflectionlist_sort_lfom(ReflectionList *in); #endif /* REFLECTION_H */ diff --git a/src/utils.c b/src/utils.c index 4618b44..1e750a8 100644 --- a/src/utils.c +++ b/src/utils.c @@ -14,6 +14,7 @@ #include <gsl/gsl_matrix.h> #include "utils.h" +#include "string.h" /* Return the MOST POSITIVE of two numbers */ unsigned int biggest(signed int a, signed int b) { @@ -81,3 +82,30 @@ void matrix_renormalise(gsl_matrix *m) { } +size_t skipspace(const char *s) { + + size_t i; + + for ( i=0; i<strlen(s); i++ ) { + if ( (s[i] != ' ') && (s[i] != '\t') ) return i; + } + + return strlen(s); + +} + +void chomp(char *s) { + + size_t i; + + if ( !s ) return; + + for ( i=0; i<strlen(s); i++ ) { + if ( (s[i] == '\n') || (s[i] == '\r') ) { + s[i] = '\0'; + return; + } + } + +} + diff --git a/src/utils.h b/src/utils.h index 1d268e6..9f8d636 100644 --- a/src/utils.h +++ b/src/utils.h @@ -18,6 +18,7 @@ #endif #include <gsl/gsl_matrix.h> +#include <math.h> extern unsigned int biggest(signed int a, signed int b); extern unsigned int smallest(signed int a, signed int b); @@ -28,6 +29,8 @@ extern double angle_between_d(double x1, double y1, double z1, double x2, double extern double lambda(double voltage); extern void matrix_renormalise(gsl_matrix *m); extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2); +extern size_t skipspace(const char *s); +extern void chomp(char *s); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) |