diff options
| -rw-r--r-- | Makefile.am | 3 | ||||
| -rw-r--r-- | doc/man/process_hkl.1 | 15 | ||||
| -rw-r--r-- | scripts/plot-cc-and-scale.R | 4 | ||||
| -rw-r--r-- | src/process_hkl.c | 94 |
4 files changed, 69 insertions, 47 deletions
diff --git a/Makefile.am b/Makefile.am index 35986f75..dbb15aa7 100644 --- a/Makefile.am +++ b/Makefile.am @@ -179,7 +179,8 @@ script_DATA = scripts/alternate-stream scripts/cell-please \ scripts/split-indexed scripts/stream_grep scripts/zone-axes \ scripts/Rsplit_surface scripts/Rsplit_surface.py \ scripts/clean-stream.py scripts/fg-graph scripts/truncate-stream \ - scripts/gen-sfs-expand scripts/add-beam-params scripts/find-pairs + scripts/gen-sfs-expand scripts/add-beam-params \ + scripts/find-pairs scripts/plot-cc-and-scale.R EXTRA_DIST += $(script_DATA) diff --git a/doc/man/process_hkl.1 b/doc/man/process_hkl.1 index b5ff885f..330ad55e 100644 --- a/doc/man/process_hkl.1 +++ b/doc/man/process_hkl.1 @@ -69,10 +69,9 @@ Stop processing after \fIn\fR crystals have been successfully merged. The defau .PD 0 .IP \fB--scale\fR .PD -Attempt to scale each pattern for best fit with the current model. +Perform a second pass through the input, scaling each crystal's intensities to best fit the initial model. -Scaling using process_hkl doesn't work very well: use \fBpartialator\fR if you -need more advanced merging techniques. +Use \fBpartialator\fR if you need more advanced merging techniques. .PD 0 .IP \fB--no-polarisation\fR @@ -106,6 +105,16 @@ Merge crystals only if they diffract to beyond \fIn\fR Angstroms resolution. Th .PD Integrate \fIn\fR nm^-1 higher than the apparent resolution limit of each individual crystal. \fIn\fR can be negative to integrate \fIlower\fR than the apparent resolution limit. The default is \fB--rescut=inf\fR, which means no resolution cutoff at all. +.PD 0 +.IP \fB--min-cc=\fIn\fR +.PD +Perform a second pass through the input, merging crystals only if their correlation with the initial model is at least \fIn\fR. + +.PD 0 +.IP \fB--stat=\fIfilename\fR +.PD +Perform a second pass through the input and, for each crystal merged, write a line to \fIfilename\fR containing the filename, scale factor and correlation coefficient with the initial model. The scale factors will all be 1 unless \fB--scale\fR is also used. + .SH CHOICE OF POINT GROUP FOR MERGING One of the main features of serial crystallography is that the orientations of diff --git a/scripts/plot-cc-and-scale.R b/scripts/plot-cc-and-scale.R index 75850e85..25df0a63 100644 --- a/scripts/plot-cc-and-scale.R +++ b/scripts/plot-cc-and-scale.R @@ -2,7 +2,7 @@ # # by Takanori Nakane (nakane.t@gmail.com) -scale_cc <- read.table("scale.csv", sep=",", col.names=c("ImageName", "ScaleFactor", "CC")) +scale_cc <- read.table("stats.dat", sep=" ", col.names=c("ImageName", "ScaleFactor", "CC")) # scatter plot # FIXME: sometimes white lines appear on the plot. @@ -20,4 +20,4 @@ abline(h=mean(sf[,2])+(-5:5)*sd(sf[,2]), col=2) hist(scale_cc$ScaleFactor, breaks=100, main="Distribution of Scale Factor", xlab="Scale Factor") abline(v=mean(sf[,2])+(-5:5)*sd(sf[,2]), col=2) -par(mfrow=c(1,1))
\ No newline at end of file +par(mfrow=c(1,1)) diff --git a/src/process_hkl.c b/src/process_hkl.c index 5bf44124..9abfc501 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -85,7 +85,6 @@ static void show_help(const char *s) " --min-snr=<n> Require individual intensity measurements to\n" " have I > n * sigma(I). Default: -infinity.\n" " --min-cc=<n> Reject frames with CC less than n. Default: infinity.\n" -" have I > n * sigma(I). Default: -infinity.\n" " --max-adu=<n> Maximum peak value. Default: infinity.\n" " --min-res=<n> Merge only crystals which diffract above <n> A.\n" " --push-res=<n> Integrate higher than apparent resolution cutoff.\n" @@ -182,10 +181,11 @@ static double scale_intensities(RefList *reference, RefList *new, return s; } + static double cc_intensities(RefList *reference, RefList *new, - const SymOpList *sym) + const SymOpList *sym) { - // x is reference + /* "x" is "reference" */ float s_xy = 0.0; float s_x = 0.0; float s_y = 0.0; @@ -234,6 +234,7 @@ static double cc_intensities(RefList *reference, RefList *new, return (s_xy - s_x*s_y/n) / sqrt(t1*t2); } + static double *check_hist_size(int n, double *hist_vals) { int ns; @@ -257,32 +258,38 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr, double **hist_vals, signed int hist_h, signed int hist_k, signed int hist_l, int *hist_n, int config_nopolar, double min_snr, double max_adu, - double push_res, double min_cc, int do_scale, FILE *stat) + double push_res, double min_cc, int do_scale, + FILE *stat) { Reflection *refl; RefListIterator *iter; + RefList *new_refl; double scale, cc; + new_refl = crystal_get_reflections(cr); + /* First, correct for polarisation */ if ( !config_nopolar ) { - polarisation_correction(crystal_get_reflections(cr), - crystal_get_cell(cr), image); + polarisation_correction(new_refl, crystal_get_cell(cr), image); } if ( reference != NULL ) { - scale = scale_intensities(reference, - crystal_get_reflections(cr), sym); - cc = cc_intensities(reference, crystal_get_reflections(cr), sym); - if (stat != NULL) { - fprintf(stat, "%s,%f,%f\n",image->filename, scale, cc); + if ( do_scale ) { + scale = scale_intensities(reference, new_refl, sym); + } else { + scale = 1.0; } - if (cc < min_cc || scale <= 0) return 1; + cc = cc_intensities(reference, new_refl, sym); + if ( cc < min_cc ) return 1; + if ( isnan(scale) ) return 1; + if ( scale <= 0.0 ) return 1; + if ( stat != NULL ) { + fprintf(stat, "%s %f %f\n", image->filename, scale, cc); + } + } else { scale = 1.0; } - if (!do_scale) scale = 1.0; - - if ( isnan(scale) ) return 1; for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; @@ -391,11 +398,12 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, int n_crystals_seen = 0; FILE *stat = NULL; - if (stat_output != NULL) { - stat = fopen(stat_output, "w"); - if (stat == NULL) { - ERROR("Failed to open statistics output file %s\n", stat_output); - } + if ( stat_output != NULL ) { + stat = fopen(stat_output, "w"); + if ( stat == NULL ) { + ERROR("Failed to open statistics output file %s\n", + stat_output); + } } do { @@ -467,8 +475,9 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, set_esd_intensity(refl, sqrt(var)/sqrt(red)); } - if (stat != NULL) - fclose(stat); + if ( stat != NULL ) { + fclose(stat); + } return 0; } @@ -504,7 +513,7 @@ int main(int argc, char *argv[]) double min_res = 0.0; double push_res = +INFINITY; double min_cc = -INFINITY; - int do_scale = 0; + int twopass = 0; /* Long options */ const struct option longopts[] = { @@ -513,7 +522,7 @@ int main(int argc, char *argv[]) {"output", 1, NULL, 'o'}, {"start-after", 1, NULL, 's'}, {"stop-after", 1, NULL, 'f'}, - {"scale", 0, &do_scale, 1}, + {"scale", 0, &config_scale, 1}, {"no-polarisation", 0, &config_nopolar, 1}, {"no-polarization", 0, &config_nopolar, 1}, {"symmetry", 1, NULL, 'y'}, @@ -549,10 +558,6 @@ int main(int argc, char *argv[]) output = strdup(optarg); break; - case 9 : - stat_output = strdup(optarg); - break; - case 's' : errno = 0; start_after = strtod(optarg, &rval); @@ -649,7 +654,12 @@ int main(int argc, char *argv[]) ERROR("Invalid value for --min-cc.\n"); return 1; } - config_scale = 1; + twopass = 1; + break; + + case 9 : + stat_output = strdup(optarg); + twopass = 1; break; case 0 : @@ -663,8 +673,6 @@ int main(int argc, char *argv[]) } - if (do_scale) config_scale = 1; - if ( filename == NULL ) { ERROR("Please specify filename using the -i option\n"); return 1; @@ -729,17 +737,21 @@ int main(int argc, char *argv[]) } + /* Need to do a second pass if we are scaling */ + if ( config_scale ) twopass = 1; + hist_i = 0; r = merge_all(st, model, NULL, sym, &hist_vals, hist_h, hist_k, hist_l, &hist_i, config_nopolar, min_measurements, min_snr, - max_adu, start_after, stop_after, min_res, push_res, min_cc, do_scale, stat_output); + max_adu, start_after, stop_after, min_res, push_res, + min_cc, config_scale, stat_output); fprintf(stderr, "\n"); if ( r ) { ERROR("Error while reading stream.\n"); return 1; } - if ( config_scale ) { + if ( twopass ) { RefList *reference; @@ -752,7 +764,7 @@ int main(int argc, char *argv[]) int r; - STATUS("Extra pass for scaling...\n"); + STATUS("Second pass for scaling and/or CCs...\n"); reference = model; model = reflist_new(); @@ -763,11 +775,12 @@ int main(int argc, char *argv[]) hist_i = 0; } - r = merge_all(st, model, reference, sym, - &hist_vals, hist_h, hist_k, hist_l, &hist_i, - config_nopolar, min_measurements, min_snr, - max_adu, start_after, stop_after, min_res, - push_res, min_cc, do_scale, stat_output); + r = merge_all(st, model, reference, sym, &hist_vals, + hist_h, hist_k, hist_l, &hist_i, + config_nopolar, min_measurements, min_snr, + max_adu, start_after, stop_after, min_res, + push_res, min_cc, config_scale, + stat_output); fprintf(stderr, "\n"); if ( r ) { ERROR("Error while reading stream.\n"); @@ -799,8 +812,7 @@ int main(int argc, char *argv[]) reflist_free(model); free(output); free(filename); - if (stat_output != NULL) - free(stat_output); + free(stat_output); return 0; } |