diff options
Diffstat (limited to 'src/render_hkl.c')
-rw-r--r-- | src/render_hkl.c | 242 |
1 files changed, 195 insertions, 47 deletions
diff --git a/src/render_hkl.c b/src/render_hkl.c index c4555242..e18ae90d 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -22,6 +22,8 @@ #include <getopt.h> #include <sys/types.h> #include <sys/wait.h> +#include <cairo.h> +#include <cairo-pdf.h> #include "utils.h" #include "reflections.h" @@ -37,76 +39,152 @@ static void show_help(const char *s) "Render intensity lists using POV-ray.\n" "\n" " -h, --help Display this help message.\n" +" --povray Render a 3D animation using POV-ray.\n" " -j <n> Run <n> instances of POV-ray in parallel.\n" +" --zone-axis Render a 2D zone axis pattern.\n" "\n"); } -int main(int argc, char *argv[]) +static void render_za(UnitCell *cell, double *ref, unsigned int *c) { - int c; - UnitCell *cell; - char *infile; - double *ref1; - signed int h, k, l; - unsigned int *c1; - FILE *fh; + cairo_surface_t *surface; + cairo_t *dctx; + double max_u, max_v, max_res, max_intensity, scale; + double sep_u, sep_v, max_r; + double as, bs, theta; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; - pid_t pids[MAX_PROC]; - float max; - int nproc = 1; - int i; + signed int h, k; + float wh, ht; - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {0, 0, NULL, 0} - }; + wh = 1024; + ht = 1024; - /* Short options */ - while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) { + surface = cairo_pdf_surface_create("za.pdf", wh, ht); - switch (c) { - case 'h' : { - show_help(argv[0]); - return 0; - } + if ( cairo_surface_status(surface) != CAIRO_STATUS_SUCCESS ) { + fprintf(stderr, "Couldn't create Cairo surface\n"); + cairo_surface_destroy(surface); + return; + } - case 'j' : { - nproc = atoi(optarg); - break; - } + dctx = cairo_create(surface); - case 0 : { - break; - } + /* Black background */ + cairo_rectangle(dctx, 0.0, 0.0, wh, ht); + cairo_set_source_rgb(dctx, 0.0, 0.0, 0.0); + cairo_fill(dctx); - default : { - return 1; - } + max_u = 0.0; max_v = 0.0; max_intensity = 0.0; + max_res = 0.0; + + /* Work out reciprocal lattice spacings and angles for this cut */ + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + theta = angle_between(asx, asy, asz, bsx, bsy, bsz); + as = modulus(asx, asy, asz) / 1e9; + bs = modulus(bsx, bsy, bsz) / 1e9; + STATUS("theta=%f\n", rad2deg(theta)); + + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + + double u, v, intensity, res; + int ct; + + ct = lookup_count(c, h, k, 0); + if ( ct < 1 ) continue; + + intensity = lookup_intensity(ref, h, k, 0) / (float)ct; + + res = resolution(cell, h, k, 0); + if ( res > max_res ) max_res = res; + + if ( intensity != 0 ) { + u = (double)h*as*sin(theta); + v = (double)h*as*cos(theta) + k*bs; + if ( fabs(u) > fabs(max_u) ) max_u = fabs(u); + if ( fabs(v) > fabs(max_v) ) max_v = fabs(v); + if ( fabs(intensity) > fabs(max_intensity) ) + max_intensity = fabs(intensity); } } + } - if ( (nproc > MAX_PROC) || (nproc < 1) ) { - ERROR("Number of processes is invalid.\n"); - return 1; + max_res /= 1e9; + max_u /= 0.5; + max_v /= 0.5; + printf("Maximum resolution is %f nm^-1\n", max_res); + + if ( max_intensity <= 0.0 ) { + max_r = 4.0; + goto out; } - infile = argv[optind]; + /* Choose whichever scaling factor gives the smallest value */ + scale = ((double)wh-50.0) / (2*max_u); + if ( ((double)ht-50.0) / (2*max_v) < scale ) { + scale = ((double)ht-50.0) / (2*max_v); + } - cell = load_cell_from_pdb("molecule.pdb"); - c1 = new_list_count(); - ref1 = read_reflections(infile, c1, NULL); - if ( ref1 == NULL ) { - ERROR("Couldn't open file '%s'\n", infile); - return 1; + sep_u = as * scale * cos(theta); + sep_v = bs * scale; + max_r = ((sep_u < sep_v)?sep_u:sep_v); + + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + + double u, v, intensity, val; + int ct; + + ct = lookup_count(c, h, k, 0); + if ( ct < 1 ) continue; + + intensity = lookup_intensity(ref, h, k, 0) / (float)ct; + val = 10.0*intensity/max_intensity; + + u = (double)h*as*sin(theta); + v = (double)h*as*cos(theta) + k*bs; + + cairo_arc(dctx, ((double)wh/2)+u*scale*2, + ((double)ht/2)+v*scale*2, max_r, 0, 2*M_PI); + + cairo_set_source_rgb(dctx, val, val, val); + cairo_fill(dctx); + + } } +out: + /* Centre marker */ + cairo_arc(dctx, (double)wh/2, + (double)ht/2, max_r, 0, 2*M_PI); + cairo_set_source_rgb(dctx, 1.0, 0.0, 0.0); + cairo_fill(dctx); + + cairo_surface_finish(surface); + cairo_destroy(dctx); +} + + +static void povray_render_animation(UnitCell *cell, double *ref, + unsigned int *c, unsigned int nproc) +{ + FILE *fh; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + pid_t pids[MAX_PROC]; + float max; + int i; + signed int h, k, l; + fh = fopen("render.pov", "w"); - fprintf(fh, "/* POV-Ray scene written by Synth3D */\n\n"); + fprintf(fh, "/* POV-Ray scene written by CrystFEL */\n\n"); fprintf(fh, "#include \"colors.inc\"\n"); fprintf(fh, "#include \"textures.inc\"\n\n"); fprintf(fh, "global_settings {\n"); @@ -232,9 +310,9 @@ int main(int argc, char *argv[]) int s; float val, p, r, g, b, trans; - if ( !lookup_count(c1, h, k, l) ) continue; + if ( !lookup_count(c, h, k, l) ) continue; - val = lookup_intensity(ref1, h, k, l); + val = lookup_intensity(ref, h, k, l); val = max-val; @@ -344,6 +422,76 @@ int main(int argc, char *argv[]) int r; waitpid(pids[i], &r, 0); } +} + + +int main(int argc, char *argv[]) +{ + int c; + UnitCell *cell; + char *infile; + double *ref; + unsigned int *cts; + int config_povray = 0; + int config_zoneaxis = 0; + unsigned int nproc = 1; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"povray", 0, &config_povray, 1}, + {"zone-axis", 0, &config_zoneaxis, 1}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hj:", longopts, NULL)) != -1) { + + switch (c) { + case 'h' : { + show_help(argv[0]); + return 0; + } + + case 'j' : { + nproc = atoi(optarg); + break; + } + + case 0 : { + break; + } + + default : { + return 1; + } + } + + } + + if ( (nproc > MAX_PROC) || (nproc < 1) ) { + ERROR("Number of processes is invalid.\n"); + return 1; + } + + infile = argv[optind]; + + cell = load_cell_from_pdb("molecule.pdb"); + cts = new_list_count(); + ref = read_reflections(infile, cts, NULL); + if ( ref == NULL ) { + ERROR("Couldn't open file '%s'\n", infile); + return 1; + } + + if ( config_povray ) { + povray_render_animation(cell, ref, cts, nproc); + } else if ( config_zoneaxis ) { + render_za(cell, ref, cts); + } else { + ERROR("Try again with either --povray or --zone-axis.\n"); + } + return 0; } |