diff options
| -rw-r--r-- | src/detector.c | 9 | ||||
| -rw-r--r-- | src/geometry-simple.tmp | 23 | ||||
| -rw-r--r-- | src/get_hkl.c | 21 | ||||
| -rw-r--r-- | src/hdf5-file.c | 28 | ||||
| -rw-r--r-- | src/image.h | 79 | ||||
| -rw-r--r-- | src/indexamajig.c | 25 | ||||
| -rw-r--r-- | src/parameters-lcls.tmp | 3 | ||||
| -rw-r--r-- | src/pattern_sim.c | 50 | ||||
| -rw-r--r-- | src/peaks.c | 63 | ||||
| -rw-r--r-- | src/sfac.c | 4 | ||||
| -rw-r--r-- | src/sfac.h | 2 |
11 files changed, 229 insertions, 78 deletions
diff --git a/src/detector.c b/src/detector.c index 78fb216f..0748b2ca 100644 --- a/src/detector.c +++ b/src/detector.c @@ -61,6 +61,7 @@ void record_image(struct image *image, int do_poisson) double total_energy, energy_density; double ph_per_e; double area; + double max_tt = 0.0; /* How many photons are scattered per electron? */ area = M_PI*pow(BEAM_RADIUS, 2.0); @@ -127,9 +128,17 @@ void record_image(struct image *image, int do_poisson) ERROR("Processed negative at %i,%i %f\n", x, y, counts); } + if ( image->twotheta[x + image->width*y] > max_tt ) { + max_tt = image->twotheta[x + image->width*y]; + } + } progress_bar(x, image->width-1, "Post-processing"); } + + STATUS("Max 2theta = %.2f deg, min d = %.2f nm (halve this to get the" + " voxel size for a synthesis)\n", + rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); } diff --git a/src/geometry-simple.tmp b/src/geometry-simple.tmp new file mode 100644 index 00000000..3472ff8a --- /dev/null +++ b/src/geometry-simple.tmp @@ -0,0 +1,23 @@ +/* Set up detector configuration */ +image.det.n_panels = 2; +image.det.panels = malloc(image.det.n_panels*sizeof(struct panel)); + +/* Upper panel */ +image.det.panels[0].min_x = 0; +image.det.panels[0].max_x = 1023; +image.det.panels[0].min_y = 512; +image.det.panels[0].max_y = 1023; +image.det.panels[0].cx = 512.0; +image.det.panels[0].cy = 502.0; +image.det.panels[0].clen = 50.0e-3; +image.det.panels[0].res = 13333.3; /* 75 micron pixel size */ + +/* Lower panel */ +image.det.panels[1].min_x = 0; +image.det.panels[1].max_x = 1023; +image.det.panels[1].min_y = 0; +image.det.panels[1].max_y = 511; +image.det.panels[1].cx = 512.0; +image.det.panels[1].cy = 522.0; +image.det.panels[1].clen = 50.0e-3; +image.det.panels[1].res = 13333.3; /* 75 micron pixel size */ diff --git a/src/get_hkl.c b/src/get_hkl.c index ae279a8e..735ae52a 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -41,7 +41,9 @@ static void show_help(const char *s) " --zone-axis Generate hk0 intensities only (and add\n" " Synth2D-style header.\n" " -i, --intensities=<file> Read intensities from file instead of\n" -" calculating them from scratch.\n" +" calculating them from scratch. You might use\n" +" this if you need to apply noise or twinning.\n" +" -p, --pdb=<file> PDB file from which to get the structure.\n" ); } @@ -114,6 +116,7 @@ int main(int argc, char *argv[]) unsigned int *cts; char *input = NULL; signed int h, k, l; + char *filename = NULL; /* Long options */ const struct option longopts[] = { @@ -124,11 +127,12 @@ int main(int argc, char *argv[]) {"twin", 0, &config_twin, 1}, {"zone-axis", 0, &config_za, 1}, {"intensities", 1, NULL, 'i'}, + {"pdb", 1, NULL, 'p'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "ht:o:i:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "ht:o:i:p:", longopts, NULL)) != -1) { switch (c) { case 'h' : { @@ -151,6 +155,11 @@ int main(int argc, char *argv[]) break; } + case 'p' : { + filename = strdup(optarg); + break; + } + case 0 : { break; } @@ -162,10 +171,14 @@ int main(int argc, char *argv[]) } - mol = load_molecule(); + if ( filename == NULL ) { + filename = strdup("molecule.pdb"); + } + + mol = load_molecule(filename); cts = new_list_count(); if ( input == NULL ) { - ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.6e-9), cts); + ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9), cts); } else { ideal_ref = read_reflections(input, cts); free(input); diff --git a/src/hdf5-file.c b/src/hdf5-file.c index 5573267b..4272cbe6 100644 --- a/src/hdf5-file.c +++ b/src/hdf5-file.c @@ -113,6 +113,10 @@ static void cleanup(hid_t fh) type = H5Iget_type(id); if ( type == H5I_GROUP ) H5Gclose(id); + if ( type == H5I_DATASET ) H5Dclose(id); + if ( type == H5I_DATATYPE ) H5Tclose(id); + if ( type == H5I_DATASPACE ) H5Sclose(id); + if ( type == H5I_ATTR ) H5Aclose(id); } } @@ -318,6 +322,11 @@ int hdf5_read(struct hdfile *f, struct image *image) { herr_t r; float *buf; + uint16_t *flags; + hid_t mask_dh; + + image->height = f->nx; + image->width = f->ny; buf = malloc(sizeof(float)*f->nx*f->ny); @@ -327,10 +336,23 @@ int hdf5_read(struct hdfile *f, struct image *image) ERROR("Couldn't read data\n"); return 1; } - image->data = buf; - image->height = f->nx; - image->width = f->ny; + + mask_dh = H5Dopen(f->fh, "/processing/hitfinder/masks", H5P_DEFAULT); + if ( mask_dh <= 0 ) { + ERROR("Couldn't open flags\n"); + } else { + flags = malloc(sizeof(uint16_t)*f->nx*f->ny); + r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL, + H5P_DEFAULT, flags); + if ( r < 0 ) { + ERROR("Couldn't read flags\n"); + image->flags = NULL; + } else { + image->flags = flags; + } + } + H5Dclose(mask_dh); /* Read wavelength from file */ image->lambda = get_wavelength(f); diff --git a/src/image.h b/src/image.h index db1b491e..70809d22 100644 --- a/src/image.h +++ b/src/image.h @@ -32,26 +32,26 @@ /* Structure describing a feature in an image */ struct imagefeature { - struct image *parent; - double x; - double y; - double intensity; + struct image *parent; + double x; + double y; + double intensity; /* Partner for this feature (in another feature list) or NULL */ - struct imagefeature_struct *partner; + struct imagefeature_struct *partner; /* Distance between this feature and its partner, if any. */ - double partner_d; + double partner_d; /* Reciprocal space coordinates (m^-1 of course) of this feature */ - double rx; - double ry; - double rz; + double rx; + double ry; + double rz; /* Internal use only */ - int valid; + int valid; - const char *name; + const char *name; }; /* An opaque type representing a list of image features */ @@ -80,46 +80,47 @@ struct reflhit { /* Structure describing an image */ struct image { - float *data; - double *twotheta; - UnitCell *indexed_cell; - UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; - int ncells; - struct detector det; - const char *filename; - struct reflhit *hits; - int n_hits; + float *data; + uint16_t *flags; + double *twotheta; + UnitCell *indexed_cell; + UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; + int ncells; + struct detector det; + const char *filename; + struct reflhit *hits; + int n_hits; - int id; + int id; - struct quaternion orientation; + struct quaternion orientation; /* Wavelength must always be given */ - double lambda; /* Wavelength in m */ + double lambda; /* Wavelength in m */ /* Incident intensity (if unknown, put 1.0) */ - double f0; - int f0_available; + double f0; + int f0_available; - int width; - int height; + int width; + int height; - ImageFeatureList *features; /* "Experimental" features */ + ImageFeatureList *features; /* "Experimental" features */ /* DirAx auto-indexing low-level stuff */ - int dirax_pty; - pid_t dirax_pid; - char *dirax_rbuffer; - int dirax_rbufpos; - int dirax_rbuflen; + int dirax_pty; + pid_t dirax_pid; + char *dirax_rbuffer; + int dirax_rbufpos; + int dirax_rbuflen; /* DirAx auto-indexing high-level stuff */ - int dirax_step; - int dirax_read_cell; - int best_acl; - int best_acl_nh; - int acls_tried[MAX_CELL_CANDIDATES]; - int n_acls_tried; + int dirax_step; + int dirax_read_cell; + int best_acl; + int best_acl_nh; + int acls_tried[MAX_CELL_CANDIDATES]; + int n_acls_tried; }; diff --git a/src/indexamajig.c b/src/indexamajig.c index 6cf9cd15..c365048a 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -260,18 +260,24 @@ static void *process_image(void *pargsv) image.hits = NULL; image.n_hits = 0; + /* View head-on (unit cell is tilted) */ + image.orientation.w = 1.0; + image.orientation.x = 0.0; + image.orientation.y = 0.0; + image.orientation.z = 0.0; + STATUS("Processing '%s'\n", image.filename); result = malloc(sizeof(*result)); if ( result == NULL ) return NULL; + result->peaks_sane = 0; + result->hit = 0; hdfile = hdfile_open(filename); if ( hdfile == NULL ) { - result->hit = 0; return result; } else if ( hdfile_set_first_image(hdfile, "/") ) { ERROR("Couldn't select path\n"); - result->hit = 0; return result; } @@ -291,17 +297,7 @@ static void *process_image(void *pargsv) if ( config_noisefilter ) { filter_noise(&image, data_for_measurement); } else { - - int x, y; - - for ( x=0; x<image.width; x++ ) { - for ( y=0; y<image.height; y++ ) { - float val; - val = image.data[x+image.width*y]; - data_for_measurement[x+image.width*y] = val; - } - } - + memcpy(data_for_measurement, image.data, data_size); } /* Perform 'fine' peak search */ @@ -329,13 +325,11 @@ static void *process_image(void *pargsv) } /* No cell at this point? Then we're done. */ - result->peaks_sane = 0; if ( image.indexed_cell == NULL ) goto done; /* Sanity check */ if ( !peak_sanity_check(&image, image.indexed_cell) ) { STATUS("Failed peak sanity check.\n"); - result->peaks_sane = 0; goto done; } else { result->peaks_sane = 1; @@ -372,6 +366,7 @@ static void *process_image(void *pargsv) done: free(image.data); + free(image.flags); free(image.det.panels); image_feature_list_free(image.features); free(image.hits); diff --git a/src/parameters-lcls.tmp b/src/parameters-lcls.tmp index 7ae2af48..1fb35065 100644 --- a/src/parameters-lcls.tmp +++ b/src/parameters-lcls.tmp @@ -12,3 +12,6 @@ /* Radius of X-ray beam */ #define BEAM_RADIUS (3.0e-6 / 2.0) + +/* Photon energy in eV */ +#define PHOTON_ENERGY (8000.0) diff --git a/src/pattern_sim.c b/src/pattern_sim.c index e09ab72a..130b6681 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -31,6 +31,7 @@ #include "peaks.h" #include "sfac.h" #include "reflections.h" +#include "parameters-lcls.tmp" static void show_help(const char *s) @@ -41,12 +42,21 @@ static void show_help(const char *s) "pulses of X-rays from a free electron laser.\n" "\n" " -h, --help Display this help message.\n" +" -p, --pdb=<file> PDB file from which to get the unit cell.\n" +" (The actual Bragg intensities come from the\n" +" intensities file)\n" " --simulation-details Show technical details of the simulation.\n" " --gpu Use the GPU to speed up the calculation.\n" "\n" " --near-bragg Output h,k,l,I near Bragg conditions.\n" " -n, --number=<N> Generate N images. Default 1.\n" " --no-images Do not output any HDF5 files.\n" +" -o, --output=<filename> Output HDF5 filename. Default: sim.h5.\n" +" If you choose to simulate more than one pattern,\n" +" the filename given will be postfixed with a\n" +" hyphen, the image number and '.h5'. In this\n" +" case, the default value is 'sim', such that the\n" +" files produced are sim-1.h5, sim-2.h5 and so on.\n" " -r, --random-orientation Use a randomly generated orientation\n" " (a new orientation will be used for each image).\n" " --powder Write a summed pattern of all images simulated by\n" @@ -166,7 +176,9 @@ int main(int argc, char *argv[]) int config_nosfac = 0; int config_gpu = 0; int config_powder = 0; + char *filename = NULL; char *grad_str = NULL; + char *outfile = NULL; GradientMethod grad; int ndone = 0; /* Number of simulations done (images or not) */ int number = 1; /* Number used for filename of image */ @@ -189,11 +201,13 @@ int main(int argc, char *argv[]) {"intensities", 1, NULL, 'i'}, {"powder", 0, &config_powder, 1}, {"gradients", 1, NULL, 'g'}, + {"pdb", 1, NULL, 'p'}, + {"output", 1, NULL, 'o'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hrn:i:g:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hrn:i:g:p:o:", longopts, NULL)) != -1) { switch (c) { case 'h' : { @@ -221,6 +235,16 @@ int main(int argc, char *argv[]) break; } + case 'p' : { + filename = strdup(optarg); + break; + } + + case 'o' : { + outfile = strdup(optarg); + break; + } + case 0 : { break; } @@ -232,6 +256,18 @@ int main(int argc, char *argv[]) } + if ( filename == NULL ) { + filename = strdup("molecule.pdb"); + } + + if ( outfile == NULL ) { + if ( n_images == 1 ) { + outfile = strdup("sim.h5"); + } else { + outfile = strdup("sim"); + } + } + if ( config_simdetails ) { show_details(); return 0; @@ -281,8 +317,8 @@ int main(int argc, char *argv[]) /* Define image parameters */ image.width = 1024; image.height = 1024; - image.lambda = ph_en_to_lambda(eV_to_J(1790.0)); /* Wavelength */ - cell = load_cell_from_pdb("molecule.pdb"); + image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */ + cell = load_cell_from_pdb(filename); image.filename = NULL; #include "geometry-lcls.tmp" @@ -372,7 +408,13 @@ int main(int argc, char *argv[]) char filename[1024]; - snprintf(filename, 1023, "results/sim-%i.h5", number); + if ( n_images != 1 ) { + snprintf(filename, 1023, "%s-%i.h5", + outfile, number); + } else { + strncpy(filename, outfile, 1023); + } + number++; /* Write the output file */ diff --git a/src/peaks.c b/src/peaks.c index 22e12949..c6b79d25 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -133,7 +133,9 @@ static void cull_peaks(struct image *image) } -static void integrate_peak(struct image *image, int xp, int yp, +/* Returns non-zero if peak has been vetoed. + * i.e. don't use result if return value is not zero. */ +static int integrate_peak(struct image *image, int xp, int yp, float *xc, float *yc, float *intensity, int do_polar) { @@ -150,6 +152,7 @@ static void integrate_peak(struct image *image, int xp, int yp, double val, sa, pix_area, Lsq, dsq, proj_area; float tt; double phi, pa, pb, pol; + uint16_t flags; /* Circular mask */ if ( x*x + y*y > lim ) continue; @@ -157,6 +160,12 @@ static void integrate_peak(struct image *image, int xp, int yp, if ( ((x+xp)>=image->width) || ((x+xp)<0) ) continue; if ( ((y+yp)>=image->height) || ((y+yp)<0) ) continue; + /* Veto this peak if we tried to integrate in a bad region */ + if ( image->flags != NULL ) { + flags = image->flags[(x+xp)+image->width*(y+yp)]; + if ( !((flags & 0x01) && (flags & 0x04)) ) return 1; + } + p = find_panel(&image->det, x+xp, y+yp); /* Area of one pixel */ @@ -203,6 +212,8 @@ static void integrate_peak(struct image *image, int xp, int yp, *yc = (float)yp; *intensity = 0; } + + return 0; } @@ -219,6 +230,7 @@ void search_peaks(struct image *image) int nrej_hot = 0; int nrej_pro = 0; int nrej_fra = 0; + int nrej_bad = 0; int nacc = 0; data = image->data; @@ -240,6 +252,7 @@ void search_peaks(struct image *image) int sx, sy; double max; unsigned int did_something; + int r; /* Overall threshold */ if ( data[x+width*y] < 800 ) continue; @@ -313,7 +326,13 @@ void search_peaks(struct image *image) /* Centroid peak and get better coordinates. * Don't bother doing polarisation correction, because the * intensity of this peak is only an estimate at this stage. */ - integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0); + r = integrate_peak(image, mask_x, mask_y, + &fx, &fy, &intensity, 0); + if ( r ) { + /* Bad region - don't detect peak */ + nrej_bad++; + continue; + } /* It is possible for the centroid to fall outside the image */ if ( (fx < 0.0) || (fx > image->width) @@ -337,8 +356,9 @@ void search_peaks(struct image *image) } } - STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame\n", - nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra); + STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame, " + "%i in bad regions.\n", + nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra, nrej_bad); cull_peaks(image); } @@ -492,6 +512,8 @@ void output_intensities(struct image *image, UnitCell *cell, int n_found; int n_indclose = 0; int n_foundclose = 0; + int n_veto = 0; + int n_veto_second = 0; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; @@ -543,17 +565,35 @@ void output_intensities(struct image *image, UnitCell *cell, &d, &idx); if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) { + int r; + /* f->intensity was measured on the filtered pattern, * so instead re-integrate using old coordinates. * This will produce further revised coordinates. */ - integrate_peak(image, f->x, f->y, &x, &y, &intensity, - !unpolar); + r = integrate_peak(image, f->x, f->y, &x, &y, + &intensity, !unpolar); + if ( r ) { + /* The original peak (which also went through + * integrate_peak(), but with the mangled + * image data) would have been rejected if it + * was in a bad region. Integration of the same + * peak included a bad region this time. */ + n_veto_second++; + continue; + } intensity = f->intensity; } else { - integrate_peak(image, hits[i].x, hits[i].y, - &x, &y, &intensity, !unpolar); + int r; + + r = integrate_peak(image, hits[i].x, hits[i].y, + &x, &y, &intensity, !unpolar); + if ( r ) { + /* Plain old ordinary peak veto */ + n_veto++; + continue; + } } @@ -592,9 +632,12 @@ void output_intensities(struct image *image, UnitCell *cell, printf("Peak statistics: %i peaks found by the peak search out of " "%i were close to indexed positions. " - "%i indexed positions out of %i were close to detected peaks\n", + "%i indexed positions out of %i were close to detected peaks.\n", n_foundclose, n_found, n_indclose, n_hits); - + printf("%i integrations using indexed locations were aborted because " + "they hit one or more bad pixels.\n", n_veto); + printf("%i integrations using peak search locations were aborted " + "because they hit one or more bad pixels.\n", n_veto_second); /* Blank line at end */ printf("\n"); @@ -311,7 +311,7 @@ static void centre_molecule(struct molecule *mol) /* Load PDB file into a memory format suitable for efficient(ish) structure * factor calculation */ -struct molecule *load_molecule() +struct molecule *load_molecule(const char *filename) { struct molecule *mol; FILE *fh; @@ -327,7 +327,7 @@ struct molecule *load_molecule() mol->reflections = NULL; mol->cell = NULL; - fh = fopen("molecule.pdb", "r"); + fh = fopen(filename, "r"); if ( fh == NULL ) { ERROR("Couldn't open PDB file\n"); return NULL; @@ -58,7 +58,7 @@ struct molecule /* This is so that the water background calculation can use it */ extern double complex get_sfac(const char *n, double s, double en); -extern struct molecule *load_molecule(void); +extern struct molecule *load_molecule(const char *filename); extern void free_molecule(struct molecule *mol); extern double *get_reflections(struct molecule *mol, double en, double res, unsigned int *counts); |