diff options
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | Makefile.in | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/Makefile.in | 5 | ||||
| -rw-r--r-- | src/cache.c | 80 | ||||
| -rw-r--r-- | src/cache.h | 34 | ||||
| -rw-r--r-- | src/control.h | 3 | ||||
| -rw-r--r-- | src/displaywindow.c | 1 | ||||
| -rw-r--r-- | src/ipr.c | 24 | ||||
| -rw-r--r-- | src/itrans.c | 3 | ||||
| -rw-r--r-- | src/main.c | 14 | ||||
| -rw-r--r-- | src/mrc.c | 4 | ||||
| -rw-r--r-- | src/peakdetect.c | 48 | ||||
| -rw-r--r-- | src/reflections.c | 8 | ||||
| -rw-r--r-- | src/reflections.h | 1 |
15 files changed, 185 insertions, 46 deletions
diff --git a/Makefile.am b/Makefile.am index 91e651a..224d9a0 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,3 +1,3 @@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src/control.h src/readpng.h src/mrc.h src/imagedisplay.h src/main.h data/displaywindow.ui src/peakdetect.h \ - src/utils.h src/itrans.h src/qdrp.h src/ipr.h + src/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h SUBDIRS = src data diff --git a/Makefile.in b/Makefile.in index 70e91b0..742d10c 100644 --- a/Makefile.in +++ b/Makefile.in @@ -159,7 +159,7 @@ target_alias = @target_alias@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src/control.h src/readpng.h src/mrc.h src/imagedisplay.h src/main.h data/displaywindow.ui src/peakdetect.h \ - src/utils.h src/itrans.h src/qdrp.h src/ipr.h + src/utils.h src/itrans.h src/qdrp.h src/ipr.h src/cache.h SUBDIRS = src data all: config.h diff --git a/src/Makefile.am b/src/Makefile.am index 1b29a7a..831d7e9 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,5 +1,5 @@ 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 peakdetect.c +dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c peakdetect.c cache.c dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" diff --git a/src/Makefile.in b/src/Makefile.in index d6994f2..cb6f316 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -47,7 +47,7 @@ am_dtr_OBJECTS = main.$(OBJEXT) displaywindow.$(OBJEXT) \ trackball.$(OBJEXT) reflections.$(OBJEXT) readpng.$(OBJEXT) \ mrc.$(OBJEXT) imagedisplay.$(OBJEXT) utils.$(OBJEXT) \ itrans.$(OBJEXT) qdrp.$(OBJEXT) ipr.$(OBJEXT) \ - peakdetect.$(OBJEXT) + peakdetect.$(OBJEXT) cache.$(OBJEXT) dtr_OBJECTS = $(am_dtr_OBJECTS) dtr_DEPENDENCIES = DEFAULT_INCLUDES = -I. -I$(top_builddir)@am__isrc@ @@ -152,7 +152,7 @@ sysconfdir = @sysconfdir@ target_alias = @target_alias@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@ -dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c peakdetect.c +dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c ipr.c peakdetect.c cache.c dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" @@ -222,6 +222,7 @@ mostlyclean-compile: distclean-compile: -rm -f *.tab.c +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cache.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/imagedisplay.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ipr.Po@am__quote@ diff --git a/src/cache.c b/src/cache.c new file mode 100644 index 0000000..1dac3ec --- /dev/null +++ b/src/cache.c @@ -0,0 +1,80 @@ +/* + * cache.c + * + * Save the reflection datablock to save having to recalculate it + * + * (c) 2007 Gordon Ball <gfb21@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#include "reflections.h" +#include "cache.h" +#include <stdlib.h> +#include <stdio.h> +#include <stdint.h> +#include <string.h> + + ReflectionContext *load_rctx_from_file(const char *filename) { + FILE *f; + CacheHeader ch; + ReflectionContext *rctx; + Reflection r; + Reflection *cr; + Reflection_NoPointer rnp; + rctx = reflection_init(); + f = fopen(filename, "rb"); + fread(&ch,sizeof(CacheHeader),1,f); + + int i; + for (i=0;i<ch.count;i++) { + + fread(&rnp,sizeof(Reflection_NoPointer),1,f); + + cr = malloc(sizeof(Reflection)); + memcpy(cr,&rnp,sizeof(Reflection)); + //printf("reading (%f,%f,%f) i=%f (%d,%d,%d) %d\n",cr->x,cr->y,cr->z,cr->intensity,cr->h,cr->k,cr->l,cr->type); + reflection_add_from_reflection(rctx,cr); + } + + fclose(f); + return rctx; + } + + int save_rctx_to_file(const char *filename, ReflectionContext *rctx) { + FILE *f; + CacheHeader ch; + Reflection *r; + Reflection_NoPointer *rnp; + char top[16] = "DTR-CACHE"; + int count=0; + + rnp = malloc(sizeof(Reflection_NoPointer)); + + r = rctx->reflections; + do { + count++; + } while ((r = r->next) != NULL ); + + + f = fopen(filename, "wb"); + memcpy(&ch.top,&top,sizeof(top)); + ch.count = count; + ch.scale = 0.; //temp, currently doesn't do anything + fwrite(&ch,sizeof(CacheHeader),1,f); + + r = rctx->reflections; + + do { + memcpy(rnp,r,sizeof(Reflection_NoPointer)); + //printf("writing (%f,%f,%f) i=%f (%d,%d,%d) %d\n",rnp->x,rnp->y,rnp->z,rnp->intensity,rnp->h,rnp->k,rnp->l,rnp->type); + fwrite(rnp,sizeof(Reflection_NoPointer),1,f); + //fwrite(r,sizeof(Reflection_NoPointer),1,f); + } while ((r = r->next) != NULL ); + + fclose(f); + return 0; + } + + +
\ No newline at end of file diff --git a/src/cache.h b/src/cache.h new file mode 100644 index 0000000..9320703 --- /dev/null +++ b/src/cache.h @@ -0,0 +1,34 @@ +#ifndef CACHE_H_ +#define CACHE_H_ + + + +typedef struct struct_cacheheader { + char top[16]; + int count; + double scale; +} CacheHeader; +/* + * Necessary to define this because I'm developing this on an amd64 machine, + * so pointers are 8 bytes instead of 4 on i386, so can't just dump reflection structs to disk +*/ +typedef struct struct_specialreflect { + double x; + double y; + double z; + double intensity; + + signed int h; + signed int k; + signed int l; + + ReflectionType type; +} Reflection_NoPointer; + +extern ReflectionContext *load_rctx_from_file(const char *filename); + +int save_rctx_to_file(const char *filename, ReflectionContext *rctx); + + +#endif /*CACHE_H_*/ + diff --git a/src/control.h b/src/control.h index d557f41..05a3c5d 100644 --- a/src/control.h +++ b/src/control.h @@ -20,7 +20,8 @@ typedef enum ift_enum { INPUT_NONE, INPUT_QDRP, - INPUT_MRC + INPUT_MRC, + INPUT_CACHE } InputFileType; typedef enum { diff --git a/src/displaywindow.c b/src/displaywindow.c index 58c9e10..8f5cd7f 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -345,6 +345,7 @@ static gint displaywindow_gl_expose(GtkWidget *widget, GdkEventExpose *event, Co reflection = ctx->reflectionctx->reflections; quadric = gluNewQuadric(); do { + if ( reflection->type == REFLECTION_CENTRAL ) { glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue); glPushMatrix(); @@ -99,6 +99,7 @@ static void ipr_try_compact(Basis *basis) { * This would probably be just as easy for the user to do... */ static Basis *ipr_choose_initial_basis(ControlContext *ctx) { + printf("ipr: c_i_b: starting\n"); Basis *basis; Reflection *reflection; int cur_max; @@ -111,10 +112,11 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) { basis->c.x = 0; basis->c.y = 0; basis->c.z = 0; basis->c.modulus = 1e30; reflection = ctx->reflectionctx->reflections->next; /* First reflection, skip initial placeholder */ cur_max = 1; + printf("ipr: c_i_b: starting loop\n"); do { double mod; unsigned int angle_works; - + //printf("ipr: c_i_b: iteration - reflection (%f,%f,%f)\n",reflection->x,reflection->y,reflection->z); mod = modulus(reflection->x, reflection->y, reflection->z); angle_works = 1; @@ -195,12 +197,16 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) { signed int h, k, l; /* NB This assumes the diffraction pattern is "vaguely" centered... */ - if ( ctx->fmode == FORMULATION_PIXELSIZE ) { - max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) * ctx->pixel_size; - } else { - max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) / (ctx->lambda * ctx->camera_length); + if ( ctx->inputfiletype != INPUT_CACHE) { + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) * ctx->pixel_size; + } else { + max_res = sqrt(ctx->width*ctx->width + ctx->height*ctx->height) / (ctx->lambda * ctx->camera_length); + } + max_res = max_res / 4; + } else { + max_res = 2e10; //works until I put some more in the reflect.cache header } - max_res = max_res / 4; max_order_a = max_res/basis->a.modulus; max_order_b = max_res/basis->b.modulus; @@ -244,20 +250,24 @@ static void ipr_sort(ReflectionContext *rctx) { } int ipr_reduce(ControlContext *ctx) { - + printf("ipr: starting\n"); Basis *basis; + printf("ipr: calling choose_initial_basis\n"); basis = ipr_choose_initial_basis(ctx); ctx->reflectionctx->basis = basis; if ( !basis ) return -1; /* Get rid of the original list and replace it with the prediction list */ + printf("ipr: calling reflection_clear\n"); reflection_clear(ctx->reflectionctx); free(ctx->reflectionctx); + printf("ipr: calling ipr_generate\n"); ctx->reflectionctx = ipr_generate(ctx, basis); ctx->reflectionctx->basis = basis; /* Sort the reflections into order of increasing g */ + printf("ipr: calling ipr_sort\n"); ipr_sort(ctx->reflectionctx); return 0; diff --git a/src/itrans.c b/src/itrans.c index c210645..27e14b6 100644 --- a/src/itrans.c +++ b/src/itrans.c @@ -515,7 +515,7 @@ static unsigned int itrans_peaksearch_iterative(int16_t *image, ControlContext * for (i=0;i<n_reflections;i++) { px = gsl_matrix_get(p,0,i); py = gsl_matrix_get(p,1,i); - printf("Reflection %d (%lf,%lf)\n",i,px,py); + //printf("Reflection %d (%lf,%lf)\n",i,px,py); if ( ctx->fmode == FORMULATION_PIXELSIZE ) { reflection_add_from_reciprocal(ctx, (px-ctx->x_centre)*ctx->pixel_size, (py-ctx->y_centre)*ctx->pixel_size, tilt_degrees, 1.0); @@ -554,5 +554,6 @@ void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degre ctx->first_image = 0; printf(" ..... %i peaks found\n", n_reflections); + } @@ -27,6 +27,7 @@ #include "mrc.h" #include "qdrp.h" #include "ipr.h" +#include "cache.h" static gint main_method_window_response(GtkWidget *method_window, gint response, ControlContext *ctx) { @@ -51,11 +52,19 @@ static gint main_method_window_response(GtkWidget *method_window, gint response, gtk_widget_destroy(method_window); while ( gtk_events_pending() ) gtk_main_iteration(); - + if ( ctx->inputfiletype == INPUT_QDRP ) { val = qdrp_read(ctx); } else if ( ctx->inputfiletype == INPUT_MRC ) { val = mrc_read(ctx); + } else if ( ctx->inputfiletype == INPUT_CACHE ) { + ctx->reflectionctx = load_rctx_from_file(ctx->filename); + val=0; + } + + if (ctx->inputfiletype == INPUT_QDRP || ctx->inputfiletype == INPUT_MRC) { + printf("Saving reflection block to reflect.cache\n"); + save_rctx_to_file("reflect.cache",ctx->reflectionctx); } if ( !val && (ctx->rmode == RECONSTRUCTION_PREDICTION) ) { @@ -161,6 +170,9 @@ int main(int argc, char *argv[]) { } else if ( strcmp(filename+(strlen(filename)-4), ".mrc") == 0 ) { printf("MRC tomography file detected.\n"); ctx->inputfiletype = INPUT_MRC; + } else if (strcmp(filename, "reflect.cache") == 0 ) { + printf("reflect.cache detected.\n"); + ctx->inputfiletype = INPUT_CACHE; } else { fprintf(stderr, "Input file must either be an MRC tomography file or a QDRP-style control file.\n"); return 1; @@ -56,6 +56,7 @@ int mrc_read(ControlContext *ctx) { /* Read all extended headers, one by one */ extsize = 4*mrc.numintegers + 4*mrc.numfloats; + printf("extsize=%d numintegers=%d numfloats=%d\n",extsize,mrc.numintegers,mrc.numfloats); if ( extsize > sizeof(MRCExtHeader) ) { fclose(fh); fprintf(stderr, "MR: MRC extended header is too big - tweak mrc.h\n"); @@ -69,7 +70,7 @@ int mrc_read(ControlContext *ctx) { printf("pixel_size=%f\n", pixel_size); nx = mrc.nx; ny = mrc.ny; - nz = (mrc.nx > mrc.ny)?mrc.nx:mrc.ny; + nz = (mrc.nx > mrc.ny)?mrc.nx:mrc.ny; //isn't this the number of images? ie usually less than nx or ny sx = nx * pixel_size; sy = ny * pixel_size; sz = nz * pixel_size; @@ -108,6 +109,7 @@ int mrc_read(ControlContext *ctx) { } } imagedisplay_mark_point(sum_id, max_x, max_y); + printf("max_x=%d max_y=%d\n",max_x,max_y); ctx->x_centre = max_x; ctx->y_centre = max_y; diff --git a/src/peakdetect.c b/src/peakdetect.c index 666e23f..37e88c5 100644 --- a/src/peakdetect.c +++ b/src/peakdetect.c @@ -22,7 +22,7 @@ void renormalise(gsl_matrix *m) { double max,min; gsl_matrix_minmax(m,&min,&max); - printf("Renormalising min=%lf max=%lf\n",min,max); + //printf("Renormalising min=%lf max=%lf\n",min,max); gsl_matrix_add_constant(m,0.-min); gsl_matrix_scale(m,1./(max-min)); } @@ -44,7 +44,7 @@ gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image) { gsl_matrix_set(raw,i,j,(double)image[i+j*ctx->width]); } } - printf("Created %dx%d matrix\n",ctx->width,ctx->height); + //printf("Created %dx%d matrix\n",ctx->width,ctx->height); renormalise(raw); return raw; } @@ -191,7 +191,7 @@ void local_sigma_filter(gsl_matrix *m, int r, double k) { int i,j; int interval = (m->size1-r)/20; //printf("lsf: starting loop\n"); - printf("lsf: "); + //printf("lsf: "); //for (i=r;i<m->size1-r;i++) { // for (j=r;j<m->size2-r;j++) { for (i=0;i<m->size1;i++) { @@ -216,9 +216,9 @@ void local_sigma_filter(gsl_matrix *m, int r, double k) { gsl_matrix_set(new,i,j,0.); } } - if (i % interval == 0) printf("."); + //if (i % interval == 0) printf("."); } - printf("done\n"); + //printf("done\n"); gsl_matrix_memcpy(m,new); gsl_matrix_free(new); } @@ -280,7 +280,7 @@ gsl_matrix* generate_flat_kernel(int half) { gsl_matrix* matrix_expand(gsl_matrix *m, int oldsize, int newsize) { gsl_matrix *new; - printf("me: %d->%d\n",oldsize,newsize); + //printf("me: %d->%d\n",oldsize,newsize); new = gsl_matrix_calloc(2,newsize); //printf("me: calloc(2,%d)\n",newsize); @@ -354,7 +354,7 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) { p = gsl_matrix_calloc(2,size); double com_x,com_y; - printf("ff: starting loop\n"); + //printf("ff: starting loop\n"); for (i=0;i<m->size1;i++) { for (j=0;j<m->size2;j++) { if (gsl_matrix_get(m,i,j)>threshold) { @@ -373,7 +373,7 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) { } com_x /= com_n; com_y /= com_n; - printf("ff: point CoM (%lf,%lf)\n",com_x,com_y); + //printf("ff: point CoM (%lf,%lf)\n",com_x,com_y); gsl_matrix_set(p,0,n,com_x); gsl_matrix_set(p,1,n,com_y); n++; @@ -385,11 +385,11 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) { } } } - printf("ff: ending loop, found %d\n",n); + //printf("ff: ending loop, found %d\n",n); *count = n; - printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); + //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); p = matrix_expand(p,size,n); - printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); + //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); return p; } @@ -408,35 +408,25 @@ gsl_matrix* iterate(gsl_matrix *m, int *count) { gsl_matrix *p; gsl_matrix *kernel; - mean = image_mean(m); - sigma = image_sigma(m,mean); - printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma); - - - printf("Iterate: generating kernel\n"); + //printf("Iterate: generating kernel\n"); kernel = generate_flat_kernel(1); printf("Iterate: performing local_sigma_filter\n"); local_sigma_filter(m,10,1.); - mean = image_mean(m); - sigma = image_sigma(m,mean); - printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma); - - - printf("Iterate: starting loop\n"); + //printf("Iterate: starting loop\n"); while (1) { - printf("Iterate: smoothing"); + //printf("Iterate: smoothing"); apply_kernel(m,kernel); - printf(" (1)"); + //printf(" (1)"); apply_kernel(m,kernel); - printf(" (2)\n"); + //printf(" (2)\n"); mean = image_mean(m); sigma = image_sigma(m,mean); - printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma); + //printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma); k = (0.5-mean)/sigma; - printf("Iterate: applying sigma_filter(%lf)\n",k); + //printf("Iterate: applying sigma_filter(%lf)\n",k); sigma_filter(m,k); - printf("Iterate: floodfilling\n"); + //printf("Iterate: floodfilling\n"); p = floodfill(m,0,&cur); printf("Iterate: %d points found\n",cur); if (old < 1.05*cur) break; diff --git a/src/reflections.c b/src/reflections.c index d2bc9db..65e88a6 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -57,7 +57,7 @@ void reflection_clear(ReflectionContext *reflectionctx) { void reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) { Reflection *new_reflection; - + //printf("Added reflection (%f,%f,%f)\n",x,y,z); new_reflection = malloc(sizeof(Reflection)); new_reflection->next = NULL; new_reflection->x = x; @@ -184,3 +184,9 @@ void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, dou reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); } + +void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) { + r->next = NULL; + rctx->last_reflection->next = r; + rctx->last_reflection = r; +} diff --git a/src/reflections.h b/src/reflections.h index 7104a27..a3a3993 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -73,5 +73,6 @@ extern void reflection_add_index(ReflectionContext *reflectionctx, signed int h, extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); extern void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); +extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r); #endif /* REFLECTION_H */ |