diff options
64 files changed, 2612 insertions, 1509 deletions
@@ -1,3 +1,36 @@ +CrystFEL version 0.4.3, 10th January 2013 +----------------------------------------- + +- compare_hkl was reworked to simplify the available figures of merit, add CC*, + and to be consistent between the overall and resolution shell values. A new + option, --intensity-shells, was also added. +- Unit cells can now be handled in any setting. +- Indexamajig will now examine the unit cell given on the command line, even if + the final cell will not be checked (i.e. if --cell-reduction=none). The + lattice type will then be given to MOSFLM, which can help a lot in finding the + right solution. +- Indexamajig now records the number of saturated reflections in the stream. +- Indexamajig now processes the peaks obtained using "--peaks=hdf5" to check for + bad detector regions, peaks too close to panel edges, etc. +- New options "--integrate-saturated", "--use-saturated" and "--integrate-found" + were added to indexamajig. +- A bug was fixed which caused systematically absent (due to centering) + reflections to be predicted and integrated. +- The colour scaling in hdfsee was improved. +- powder_plot was moved to the extra programs repository. +- The configure script now checks for HDF5. +- Forbidden (due to centering) reflections are now taken into account by + check_hkl when calculating the completeness. +- The speed of pattern processing was increased significantly (in many cases) by + avoiding an unnecessary iteration. +- A crystfel.pc file was added, allowing use of "pkg-config --cflags crystfel" + and so on. +- The new option "--min-measurements" was added to process_hkl. +- The wavelength can now be specified in the beam file as an HDF5 path or as an + explicit number, just like the camera length. There is no longer any + "fallback" to a nominal value (Chun Hong Yoon). + + CrystFEL version 0.4.2, 2nd October 2012 ---------------------------------------- diff --git a/Makefile.am b/Makefile.am index bf353e2e..422a34fb 100644 --- a/Makefile.am +++ b/Makefile.am @@ -4,7 +4,7 @@ SUBDIRS = lib doc/reference libcrystfel ACLOCAL_AMFLAGS = -I m4 bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ - src/compare_hkl src/partialator src/check_hkl src/partial_sim + src/compare_hkl src/partialator src/check_hkl src/partial_sim noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_gradient_check tests/symmetry_check \ @@ -61,10 +61,6 @@ In addition, there is also: - get_hkl, for doing various simple operations on reflection lists. - - powder_plot, for turning image or reflection data into a histogram of - reciprocal space vector moduli and intensities (i.e. a - 1D powder diffraction trace). - - compare_hkl, for working out the differences (e.g. a q-dependent scaling factor) between two lists of reflections. diff --git a/configure.ac b/configure.ac index a148942a..a3a15482 100644 --- a/configure.ac +++ b/configure.ac @@ -1,6 +1,6 @@ dnl Process this file with autoconf to produce a configure script. -AC_INIT(crystfel, 0.4.2, taw@physics.org) +AC_INIT(crystfel, 0.4.3, taw@physics.org) AC_CONFIG_AUX_DIR([build-aux]) AM_CONFIG_HEADER(config.h) AM_INIT_AUTOMAKE([subdir-objects]) diff --git a/doc/examples/lcls-cxi-9keV.beam b/doc/examples/lcls-cxi-9keV.beam index d0476cd3..5241e04f 100644 --- a/doc/examples/lcls-cxi-9keV.beam +++ b/doc/examples/lcls-cxi-9keV.beam @@ -4,8 +4,8 @@ beam/fluence = 1.0e12 ; Radius of X-ray beam beam/radius = 10.0e-6 -; Photon energy in eV -beam/photon_energy = 9340.0 +; Photon energy in eV (9340 eV, but use individual value from file) +beam/photon_energy = /LCLS/photon_energy_eV ; Bandwidth: FWHM(wavelength) over wavelength. ; Note: current simulation code just uses a rectangular diff --git a/doc/examples/lcls-dec.beam b/doc/examples/lcls-dec.beam index 1f9fe3f6..358f8d44 100644 --- a/doc/examples/lcls-dec.beam +++ b/doc/examples/lcls-dec.beam @@ -4,8 +4,8 @@ beam/fluence = 2.0e11 ; Radius of X-ray beam beam/radius = 1.5e-6 -; Photon energy in eV -beam/photon_energy = 2000.0 +; Photon energy in eV (2 keV) +beam/photon_energy = /LCLS/photon_energy_eV ; Bandwidth: FWHM(wavelength) over wavelength. ; Note: current simulation code just uses a rectangular diff --git a/doc/examples/lcls-june.beam b/doc/examples/lcls-june.beam index 7d95c820..57d5e802 100644 --- a/doc/examples/lcls-june.beam +++ b/doc/examples/lcls-june.beam @@ -5,7 +5,7 @@ beam/fluence = 8.0e12 beam/radius = 1.5e-6 ; Photon energy in eV -beam/photon_energy = 2000.0 +beam/photon_energy = /LCLS/photon_energy_eV ; Bandwidth: FWHM(wavelength) over wavelength. ; Note: current simulation code just uses a rectangular diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1 index bddbc66c..889cd9cd 100644 --- a/doc/man/compare_hkl.1 +++ b/doc/man/compare_hkl.1 @@ -98,6 +98,13 @@ Discard reflections with lower than this value of 1/d in m^-1. .PD Discard reflections with higher than this value of 1/d in m^-1. +.PD 0 +.IP \fB--intensity-shells\fR +.PD +Use intensity shells instead of resolution shells. The range of shells will start at the intensity of the least intense reflection, and extend from there by 1/5000th of the intensity of the difference between the strongest and weakest reflection in the first reflection list. The pairs of reflections will also be assigned their bins according to the intensity of the reflection in the first list. +.sp +Because of the hardcoded factor of 1/5000, needed to avoid a very uneven distribution of the number of reflection pairs in each bin, you are advised not to draw strong conclusions from the results of using this option. + .SH AUTHOR This page was written by Thomas White. diff --git a/doc/man/crystfel.7 b/doc/man/crystfel.7 index 3efb3004..b2876b59 100644 --- a/doc/man/crystfel.7 +++ b/doc/man/crystfel.7 @@ -44,9 +44,6 @@ In addition, there is also: .IP \fBget_hkl\fR A tool for manipulating reflection lists, such as performing symmetry expansion. -.IP \fBpowder_plot\fR -A tool for the calculation of one-dimensional "powder" traces. - .IP \fBcompare_hkl\fR and \fBcheck_hkl\fR Tools for calculating figures of merit, such as completeness and R-factors. @@ -171,7 +168,6 @@ You should have received a copy of the GNU General Public License along with Cry .BR compare_hkl (1), .BR check_hkl (1), .BR render_hkl (1), -.BR powder_plot (1), .BR hdfsee (1), .BR get_hkl (1), .BR crystfel_geometry (5). diff --git a/doc/man/crystfel_geometry.5 b/doc/man/crystfel_geometry.5 index ccc6fad1..216aec15 100644 --- a/doc/man/crystfel_geometry.5 +++ b/doc/man/crystfel_geometry.5 @@ -15,7 +15,7 @@ See below for information about CrystFEL's beam description files. .SH CRYSTFEL DETECTOR GEOMETRY FILES The detector geometry is taken from a text file rather than hardcoded into the program. Programs which care about the geometry (particularly -\fBindexamajig\fR, \fBpattern_sim\fR and \fBpowder_plot\fR) take an argument +\fBindexamajig\fR and \fBpattern_sim\fR take an argument \fB--geometry=\fR\fIfilename\fR, where \fIfilename\fR contains the geometry. .PP A flexible (and pedantic) representation of the detector has been developed to diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1 index 3519dc5f..11da2445 100644 --- a/doc/man/indexamajig.1 +++ b/doc/man/indexamajig.1 @@ -1,7 +1,8 @@ .\" .\" indexamajig man page .\" -.\" Copyright © 2012 Thomas White <taw@physics.org> +.\" Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. .\" .\" Part of CrystFEL - crystallography with a FEL .\" @@ -12,25 +13,23 @@ indexamajig \- bulk indexing and data reduction program .SH SYNOPSIS .PP .BR indexamajig -\fB-i\fR \fIfilename\fR \fB-o\fR \fIoutput.stream\fR \fB-g\fR \fIdetector.geom\fR \fB-b\fR \fIbeamline.beam\fR \fB--peaks=\fR\fImethod\fR \fB--indexing=\fR\fImethod\fR \fB--cell-reduction=\fR\fImethod\fR +\fB-i\fR \fIfilename\fR \fB-o\fR \fIoutput.stream\fR \fB-g\fR \fIdetector.geom\fR \fB-b\fR \fIbeamline.beam\fR \fB--peaks=\fR\fImethod\fR \fB--indexing=\fR\fImethod\fR [\fBoptions\fR] \fB...\fR .PP \fBindexamajig --help\fR .SH DESCRIPTION -indexamajig takes as input a list of diffraction image files, currently in HDF5 format. For each image, it attempts to find peaks and then index the pattern. If successful, it will measure the intensities of the peaks at Bragg locations and produce a list of integrated peaks and intensities (and so on) in an output text file known as a "stream". +\fBindexamajig\fR takes a list of diffraction snapshots from crystals in random orientations and attempts to find peaks, index and integrate each one. The input is a list of diffraction image files in HDF5 format and some auxiliary files and parameters. The output is a long text file ('stream') containing the results from each image in turn. -For minimal basic use, you need to provide the list of diffraction patterns, the method which will be used to index, a file describing the geometry of the detector, a PDB file which contains the unit cell which will be used for the indexing, and that you'd like the program to output a list of intensities for each successfully indexed pattern. Here is what the minimal use might look like on the command line: +For minimal basic use, you need to provide the list of diffraction patterns, the method which will be used to index, a file describing the geometry of the detector, and a PDB file which contains the unit cell which will be used for the indexing. Here is what the minimal use might look like on the command line: .IP \fBindexamajig\fR .PD --i mypatterns.lst -j 10 -g mygeometry.geom --indexing=mosflm,dirax --peaks=hdf5 --cell-reduction=reduce -b myxfel.beam -o test.stream -p mycell.pdb --record=integrated +-i mypatterns.lst -j 10 -g mygeometry.geom --indexing=mosflm,dirax --peaks=hdf5 -b myxfel.beam -o test.stream -p mycell.pdb .PP -More typical use includes all the above, but might also include a noise or -common mode filter (--filter-noise or --filter-cm respectively) if detector -noise causes problems for the peak detection. The HDF5 files might be in some +More typical use includes all the above, but might also include extra parameters to modify the behaviour. The HDF5 files might be in some folder a long way from the current directory, so you might want to specify a full pathname to be added in front of each filename. You'll probably want to run more than one indexing job at a time (-j <n>). @@ -43,97 +42,76 @@ array, where the first two columns contain fast scan and slow scan coordinates at the given location. The value will be spread in a small cross centred on that location. -See \fBman crystfel_geometry\fR for information about how to create a geometry description file. - -You can control what information is included in the output stream using -\fB--record\fR=\fIflags\fR. Possible flags are: - -.IP \fBintegrated\fR -.PD -Include a list of reflection intensities, produced by integrating around predicted peak locations. - -.IP \fBpeaks\fR -.PD -Include peak locations and intensities from the peak search. - -.IP \fBpeaksifindexed\fR -.PD -As 'peaks', but the peak information will be written only if the pattern could be indexed. - -.IP \fBpeaksifnotindexed\fR -.PD -As 'peaks', but the peak information will be written only if the pattern could \fInot\fR be indexed. - -.PP -The default is \fB--record=integrated\fR. - -.PP -If you just want the integrated intensities of indexed peaks, use \fB--record=integrated\fR. If you just want to check that the peak detection is working, used \fB--record=peaks\fR. If you want the integrated peaks for the indexable patterns, but also want to check the peak detection for the patterns -which could not be indexed, you might use \fB--record=integrated,peaksifnotindexed\fR and then use \fBcheck-peak-detection\fR from the scripts folder to visualise the results of the peak detection. +See \fBman crystfel_geometry\fR for information about how to create a geometry description file and a beam parameters file. .SH PEAK DETECTION You can control the peak detection on the command line. Firstly, you can choose the peak detection method using \fB--peaks=\fR\fImethod\fR. Currently, two values for "method" are available. \fB--peaks=hdf5\fR will take the peak locations from the HDF5 file. It expects a two dimensional array at where size in the first dimension is the number of peaks and the size in the second dimension is three. The first two columns contain the fast scan and slow scan coordinates, the third contains the intensity. However, the intensity will be ignored since the pattern will always be re-integrated using the unit cell provided by the indexer on the basis of the peaks. You can tell indexamajig where to find this table inside each HDF5 file using \fB--hdf5-peaks=\fR\fIpath\fR. -If you use \fB--peaks=zaef\fR, indexamajig will use a simple gradient search after Zaefferer (2000). You can control the overall threshold and minimum gradient for finding a peak using \fB--threshold\fR and \fB--min-gradient\fR. Both of these have arbitrary units matching the pixel values in the data. +If you use \fB--peaks=zaef\fR, indexamajig will use a simple gradient search after Zaefferer (2000). You can control the overall threshold and minimum squared gradient for finding a peak using \fB--threshold\fR and \fB--min-gradient\fR. The threshold has arbitrary units matching the pixel values in the data, and the minimum gradient has the equivalent squared units. -A minimum peak separation can also be provided in the geometry description file -(see \fBman crystfel_geometry\fR for details). This number serves two purposes. Firstly, it is the maximum distance allowed between the peak summit and the foot point (where the gradient exceeds the minimum gradient). Secondly, it is the minimum distance allowed between one peak and another, before the later peak will be rejected. +Peaks will be rejected if the 'foot point' is further away from the 'summit' of the peak by more than the inner integration radius (see below). They will also be rejected if the peak is closer than twice the inner integration radius from another peak. You can suppress peak detection altogether for a panel in the geometry file by specifying the "no_index" value for the panel as non-zero. .SH INDEXING METHODS -You can choose between a variety of indexing methods. You can choose more than one method, in which case each method will be tried in turn until the later cell reduction step says that the cell is a "hit". Choose from: - -.IP \fBnone\fR -.PD -No indexing, peak detection only. +You can choose between a variety of indexing methods. You can choose more than one method, in which case each method will be tried in turn until one of them reports that the pattern has been successfully indexed. Choose from: .IP \fBdirax\fR .PD -DirAx will be executed. DirAx must be installed and in your PATH for this to work. Test by running \fBdirax\fR on the command line immediately before running \fBindexamajig\fR - you should see a welcome message. If not, refer to the installation instructions for DirAx. +Invoke DirAx, check linear combinations of the resulting cell axes for agreement with your cell, and then check that the cell accounts for at least half of the peaks from the peak search. .IP \fBmosflm\fR .PD -MOSFLM will be executed. MOSFLM must be installed and in your PATH for this to work. Test by running \fBipmosflm\fR on the command line immediately before running \fBindexamajig\fR - you should see a welcome message. If not, refer to the installation instructions for CCP4. +As \fBdirax\fR, but invoke MOSFLM instead. If you provide a PDB file (with \fB-p\fR), the lattice type and centering information will be passed to MOSFLM, which will then return solutions which match. Note that the lattice parameter information will \fBnot\fR be given to MOSFLM, because it has no way to make use of it. .IP \fBreax\fR .PD -An \fIexperimental\fR DPS-style algorithm will be used, which searches for a lattice with cell parameters close to those contained in the CRYST1 line of the PDB file given with \fB-p\fR or \fB--pdb\fR. If you use this option, you should also use \fB--cell-reduction=none\fR. +Run the DPS algorithm, looking for the axes of your cell. .PP -The default is \fB--indexing=none\fR. +You can add one or more of the following to the above indexing methods: +.IP \fB-raw\fR +.PD +Do not check the resulting unit cell. This option is useful when you need to determine the unit cell ab initio. Use with 'dirax' and 'mosflm' - the other indexing methods need the unit cell as input in any case, and cannot determine the unit cell ab initio. See \fB-comb\fR and \fB-axes\fR. + +.IP \fB-axes\fR +.PD +Check permutations of the axes for correspondence with your cell, but do not check linear combinations. This is useful to avoid a potential problem when one of the unit cell axis lengths is close to a multiple of one of the others. Can be used with \fBdirax\fR and \fBmosflm\fR. See \fB-raw\fR and \fB-comb\fR. -.SH CELL REDUCTION -Neither MOSFLM nor DirAx are currently able to simply search for the orientation of a known unit cell. Instead, they determine the unit cell parameters ab initio. To check if the parameters match the known unit cell, which you provide with \fB-p\fR \fIunitcell.pdb\fR or \fB--pdb=\fR\fIunitcell.pdb\fR. There is a choice of ways in which this comparison can be performed, which you can choose between using \fB--cell-reduction=\fR\fImethod\fR. The choices for \fImethod\fR are: +.IP \fB-comb\fR +.PD +Check linear combinations of the unit cell basis vectors to see if a cell can be produced which looks like your unit cell. This is the default behaviour for \fBdirax\fR and \fBmosflm\fR. See \fB-raw\fR and \fB-axes\fR. -.IP \fBnone\fR +.IP \fB-bad\fR .PD -The raw cell from the autoindexer will be used. The cell probably won't match the target cell, but it'll still get used. Use this option to test whether the patterns are basically "indexable" or not, or if you don't know the cell parameters. In the latter case, you'll need to plot a histogram of the resulting parameters from the output stream to see which are the most popular. +Do not check that the cell accounts for any of the peaks as described in \fBdirax\fR above. Might be useful to debug initial indexing problems, or if there are many multi-crystal patterns and the indexing method has no concept of multiple crystals per pattern (which, at the moment, means all of them). Can be used with any indexing method, but is generally a bad idea. -.IP \fBreduce\fR +.IP \fB-nolatt\fR .PD -Linear combinations of the raw cell will be checked against the target cell. If at least one candidate is found for each axis of the target cell, the angles will be checked to correspondence. If a match is found, this cell will be used for further processing. This option should generate the most matches. +Do not use the lattice type information from the PDB file to help guide the indexing. Use with \fBmosflm\fR, which is the only indexing method which can optionally take advantage of this information. This is the default behaviour for \fBdirax\fR. This option makes no sense for \fBreax\fR, which is intrinsically based on using known lattice information. -.IP \fBcompare\fR +.IP \fB-latt\fR .PD -The cell will be checked for correspondence after only a simple permutation of the axes. This is useful when the target cell is subject to reticular twinning, such as if one cell axis length is close to twice another. With \fB--cell-reduction=reduce\fR there would be a possibility that the axes might be confused in this situation. This happens for lysozyme (1VDS), so watch out. +This is the opposite of \fB-nolatt\fR, and is the default behaviour for \fBmosflm\fR and \fBreax\fR. .PP -The tolerance for matching with \fBreduce\fR and \fBcompare\fR can be set using \fB--tolerance=\fR\fIa,b,c,angl\fR \- see below for more details. Cells from these reduction routines are further constrained to be right-handed, but no such constraint will be applied if you use \fB--cell-reduction=none\fR. Always using a right-handed cell means that the Bijvoet pairs can be told apart. +The default indexing method is 'none', which means no indexing will be done. This is useful if you just want to check that the peak detection is working properly. .PP -If the unit cell is centered (i.e. if the space group begins with I, F, H, A, B, C), you should be careful when using "compare" for the cell reduction, since (for example) DirAx will always find a primitive unit cell, and this cell must be converted to the non-primitive conventional cell from the PDB. +Your indexing methods will be checked for validity, incompatible flags removed, and warnings given about duplicates. For example, \fBmosflm\fR and \fBmosflm-comb-latt\fR represent the same indexing method, because \fB-comb\fR and \fB-latt\fR are the default behaviour for \fBmosflm\fR. The 'long version' of each of your indexing methods will be listed in the output, and the stream will contain a record of which indexing method successfully indexed each pattern. .PP -The default is \fB--cell-reduction=none\fR. +It's risky to use \fBmosflm-nolatt\fR in conjunction with either \fB-comb\fR or \fB-axes\fR when you have a rhombohedral cell. This would be an odd thing to do anyway: why withhold the lattice information from MOSFLM if you know what it is, and want to use it to check the result? It's risky because MOSFLM will by default return the "H centered" lattice for your rhombohedral cell, and it's not completely certain that MOSFLM consistently uses one or other of the two possible conventions for the relationship between the "H" and "R" cells. It is, however, very likely that it does. + +Examples of indexing methods: 'dirax,mosflm,reax', 'dirax-raw,mosflm-raw', 'dirax-raw-bad'. + .SH PEAK INTEGRATION -If the pattern could be successfully indexed and the cell reduction gave a -positive match, peaks will be predicted in the pattern and their intensities +If the pattern could be successfully indexed, peaks will be predicted in the pattern and their intensities measured. The peak integration relies on knowing an accurate radius to integrate the peak within, and that the annulus used to estimate the background level is not so big that that it covers neighbouring peaks. However, @@ -179,11 +157,6 @@ Find peaks in the images using \fImethod\fR. See the second titled \fBPEAK DETE Index the patterns using \fImethod\fR. See the section titled \fBINDEXING METHODS\fR (above) for more information. .PD 0 -.IP \fB--cell-reduction=\fR\fImethod\fR -.PD -Use \fImethod\fR for unit cell reduction. See the section titled \fBCELL REDUCTION\fR (above) for more information. - -.PD 0 .IP "\fB-g\fR \fIfilename\fR" .IP \fB--geometry=\fR\fIfilename\fR .PD @@ -235,7 +208,7 @@ When using \fB--peaks=hdf5\fR, read the peak locations from a table in the HDF5 .PD 0 .IP \fB--tolerance=\fR\fItol\fR .PD -Set the tolerances for unit cell comparison. \fItol\fR takes the form \fIa\fR,\fIb\fR,\fIc\fR,\fIang\fR. \fIa\fR, \fIb\fR and \fIc\fR are the tolerances, in percent, for the respective direct space axes when using \fBcompare\fR or \fBcompare_ab\fR for cell reduction (see above). \fIang\fR is the tolerance in degrees for the angles. They represent the respective \fIreciprocal\fR space parameters when using \fB--cell-reduction=reduce\fR. +Set the tolerances for unit cell comparison. \fItol\fR takes the form \fIa\fR,\fIb\fR,\fIc\fR,\fIang\fR. \fIa\fR, \fIb\fR and \fIc\fR are the tolerances, in percent, for the respective direct space axes when using \fB-axes\fR in the indexing method (see below). \fIang\fR is the tolerance in degrees for the angles. When \fBnot\fR using \fB-axes\fR, they represent the respective \fIreciprocal\fR space parameters. Sorry for the horrifying inconsistency. .PD The default is \fB--tolerance=5,5,5,1.5\fR. @@ -309,11 +282,6 @@ with old scripts because this option selects the behaviour which is now the default. .PD 0 -.IP \fB--insane\fR -.PD -Normally, indexamajig will check that the projected reciprocal space positions of peaks found by the peak detection are close to reciprocal lattice points. This option disables this check, and normally is \fInot\fR a good idea. - -.PD 0 .IP \fB--no-bg-sub\fR .PD Don't subtract local background estimates from integrated intensities. This is almost never useful, but might help to detect problems from troublesome background noise. @@ -334,6 +302,11 @@ Normally, reflections which contain one or more pixels above max_adu (defined in .PD When using \fB--peaks=hdf5\fR, the peaks will be put through the same checks as if you were using \fB--peaks=zaef\fR. These checks reject peaks which are too close to panel edges, are saturated (unless you use \fB--use-saturated\fR), fall short of the minimum SNR value given by \fB--min-snr\fR, have other nearby peaks (closer than twice the inner integration radius, see \fB--int-radius\fR), or have any part in a bad region. Using this option skips this validation step, and uses the peaks directly. +.PD 0 +.IP \fB--integrate-found\fR +.PD +Without this option, peaks will be predicted in each pattern based on the crystal orientation from autoindexing and the parameters in the beam description file. With this option, indices will be assigned to the peaks found by the peak search, and the positions of those peaks used for the final integration stage. + .SH BUGS Don't run more than one indexamajig jobs simultaneously in the same working directory - they'll overwrite each other's DirAx or MOSFLM files, causing subtle @@ -348,7 +321,7 @@ This page was written by Thomas White. Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>. .SH COPYRIGHT AND DISCLAIMER -Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. +Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. .P indexamajig, and this manual, are part of CrystFEL. .P diff --git a/doc/man/process_hkl.1 b/doc/man/process_hkl.1 index 1b07dc04..1a8bc05b 100644 --- a/doc/man/process_hkl.1 +++ b/doc/man/process_hkl.1 @@ -90,8 +90,7 @@ need more advanced merging techniques. .PD Sum the intensities, instead of averaging them. This might be useful for comparing results to radially summed powder patterns, but usually gives less -useful final intensities. \fBpowder_plot\fR has more advanced features for -calculating 'powder' patterns. +useful final intensities. .PD 0 .IP \fB--max-only\fR @@ -142,5 +141,4 @@ You should have received a copy of the GNU General Public License along with Cry .BR compare_hkl (1), .BR check_hkl (1), .BR render_hkl (1), -.BR powder_plot (1), .BR partialator (1) diff --git a/doc/reference/CrystFEL-docs.sgml b/doc/reference/CrystFEL-docs.sgml index 348169bd..4e6cc877 100644 --- a/doc/reference/CrystFEL-docs.sgml +++ b/doc/reference/CrystFEL-docs.sgml @@ -8,7 +8,7 @@ <bookinfo> <title>CrystFEL Reference Manual</title> <releaseinfo> - For libcrystfel from CrystFEL 0.4.2. + For libcrystfel from CrystFEL 0.4.3. </releaseinfo> <abstract> This is the internal documentation for CrystFEL. Unless you are looking at @@ -38,6 +38,11 @@ </chapter> <chapter> + <title>Crystals</title> + <xi:include href="xml/crystal.xml"/> + </chapter> + + <chapter> <title>Statistics and R-factors</title> <xi:include href="xml/statistics.xml"><xi:fallback /></xi:include> </chapter> diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt index 18a20eb9..5df1f2d5 100644 --- a/doc/reference/CrystFEL-sections.txt +++ b/doc/reference/CrystFEL-sections.txt @@ -287,6 +287,13 @@ largest_q </SECTION> <SECTION> +<FILE>crystal</FILE> +Crystal +crystal_new +crystal_free +</SECTION> + +<SECTION> <FILE>hdf5-file</FILE> hdf5_read hdf5_write diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 22009223..ae315594 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -8,7 +8,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ - src/cell-utils.c src/integer_matrix.c \ + src/cell-utils.c src/integer_matrix.c src/crystal.c \ src/grainspotter.c if HAVE_FFTW @@ -25,7 +25,8 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ src/render.h src/index.h src/image.h \ src/filters.h src/dirax.h src/mosflm.h \ src/index-priv.h src/reax.h src/cell-utils.h \ - src/integer_matrix.h src/grainspotter.h + src/integer_matrix.h src/crystal.h \ + src/grainspotter.h AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall AM_CPPFLAGS += -I$(top_srcdir)/lib @LIBCRYSTFEL_CFLAGS@ diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h index 8212811b..de777deb 100644 --- a/libcrystfel/src/beam-parameters.h +++ b/libcrystfel/src/beam-parameters.h @@ -3,11 +3,11 @@ * * Beam parameters * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010,2012 Thomas White <taw@physics.org> + * 2010,2012-2013 Thomas White <taw@physics.org> * 2012 Chunhong Yoon * * This file is part of CrystFEL. @@ -34,6 +34,8 @@ #include <config.h> #endif +struct beam_params; + #include "hdf5-file.h" struct beam_params diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 3ec9a9d6..ee51622b 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -424,14 +424,18 @@ UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t) tt = uncentering_transformation(in, &new_centering, &new_latt); if ( tt == NULL ) return NULL; - if ( t != NULL ) *t = tt; - out = cell_transform(in, tt); if ( out == NULL ) return NULL; cell_set_lattice_type(out, new_latt); cell_set_centering(out, new_centering); + if ( t != NULL ) { + *t = tt; + } else { + tfn_free(tt); + } + return out; } diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c new file mode 100644 index 00000000..664d933e --- /dev/null +++ b/libcrystfel/src/crystal.c @@ -0,0 +1,224 @@ +/* + * crystal.c + * + * A class representing a single crystal + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2013 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "crystal.h" +#include "utils.h" + + +/** + * SECTION:crystal + * @short_description: Crystal + * @title: Crystal + * @section_id: + * @see_also: + * @include: "crystal.h" + * @Image: + * + * This structure represents a single crystal. + */ + + +struct _crystal +{ + /* The image containing the crystal */ + struct image *image; + + /* Information about the crystal */ + UnitCell *cell; + double m; /* Mosaicity in radians */ + double osf; + double profile_radius; + int pr_dud; + double resolution_limit; + + /* Integrated (or about-to-be-integrated) reflections */ + RefList *reflections; + long long int n_saturated; /* Number of overloads */ + + /* User flag, e.g. for "this is a bad crystal". */ + int user_flag; +}; + + +/************************** Constructor / Destructor **************************/ + + +/** + * crystal_new: + * + * Create a new %Crystal. + * + * Returns: the new unit cell, or NULL on failure. + * + */ +Crystal *crystal_new() +{ + Crystal *cryst; + + cryst = malloc(sizeof(Crystal)); + if ( cryst == NULL ) return NULL; + + cryst->cell = NULL; + cryst->reflections = NULL; + cryst->resolution_limit = 0.0; + cryst->n_saturated = 0; + + return cryst; +} + + +/** + * crystal_free: + * @cryst: A %Crystal to free. + * + * Frees a %Crystal, and all internal resources concerning that crystal. + * + */ +void crystal_free(Crystal *cryst) +{ + if ( cryst == NULL ) return; + free(cryst); +} + + +/********************************** Getters ***********************************/ + + +UnitCell *crystal_get_cell(Crystal *cryst) +{ + return cryst->cell; +} + + +double crystal_get_profile_radius(Crystal *cryst) +{ + return cryst->profile_radius; +} + + +RefList *crystal_get_reflections(Crystal *cryst) +{ + return cryst->reflections; +} + + +double crystal_get_resolution_limit(Crystal *cryst) +{ + return cryst->resolution_limit; +} + + +long long int crystal_get_num_saturated_reflections(Crystal *cryst) +{ + return cryst->n_saturated; +} + + +struct image *crystal_get_image(Crystal *cryst) +{ + return cryst->image; +} + + +double crystal_get_osf(Crystal *cryst) +{ + return cryst->osf; +} + + +int crystal_get_user_flag(Crystal *cryst) +{ + return cryst->user_flag; +} + + +double crystal_get_mosaicity(Crystal *cryst) +{ + return cryst->m; +} + + +/********************************** Setters ***********************************/ + + +void crystal_set_cell(Crystal *cryst, UnitCell *cell) +{ + cryst->cell = cell; +} + + +void crystal_set_profile_radius(Crystal *cryst, double r) +{ + cryst->profile_radius = r; +} + + +void crystal_set_reflections(Crystal *cryst, RefList *reflist) +{ + cryst->reflections = reflist; +} + + +void crystal_set_resolution_limit(Crystal *cryst, double res) +{ + cryst->resolution_limit = res; +} + + +void crystal_set_num_saturated_reflections(Crystal *cryst, long long int n) +{ + cryst->n_saturated = n; +} + + +void crystal_set_image(Crystal *cryst, struct image *image) +{ + cryst->image = image; +} + + +void crystal_set_osf(Crystal *cryst, double osf) +{ + cryst->osf = osf; +} + + +void crystal_set_user_flag(Crystal *cryst, int user_flag) +{ + cryst->user_flag = user_flag; +} + + +void crystal_set_mosaicity(Crystal *cryst, double m) +{ + cryst->m = m; +} diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h new file mode 100644 index 00000000..3d0ad9d1 --- /dev/null +++ b/libcrystfel/src/crystal.h @@ -0,0 +1,74 @@ +/* + * crystal.h + * + * A class representing a single crystal + * + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2013 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef CRYSTAL_H +#define CRYSTAL_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include "cell.h" +#include "reflist.h" + + +/** + * Crystal: + * + * This data structure is opaque. You must use the available accessor functions + * to read and write its contents. + **/ +typedef struct _crystal Crystal; + + +extern Crystal *crystal_new(void); +extern void crystal_free(Crystal *cryst); + +extern UnitCell *crystal_get_cell(Crystal *cryst); +extern double crystal_get_profile_radius(Crystal *cryst); +extern RefList *crystal_get_reflections(Crystal *cryst); +extern double crystal_get_resolution_limit(Crystal *cryst); +extern long long int crystal_get_num_saturated_reflections(Crystal *cryst); +extern int crystal_get_user_flag(Crystal *cryst); +extern double crystal_get_osf(Crystal *cryst); +extern struct image *crystal_get_image(Crystal *cryst); +extern double crystal_get_mosaicity(Crystal *cryst); + +extern void crystal_set_cell(Crystal *cryst, UnitCell *cell); +extern void crystal_set_profile_radius(Crystal *cryst, double r); +extern void crystal_set_reflections(Crystal *cryst, RefList *reflist); +extern void crystal_set_resolution_limit(Crystal *cryst, double res); +extern void crystal_set_num_saturated_reflections(Crystal *cryst, + long long int n); +extern void crystal_set_user_flag(Crystal *cryst, int flag); +extern void crystal_set_osf(Crystal *cryst, double osf); +extern void crystal_set_image(Crystal *cryst, struct image *image); +extern void crystal_set_mosaicity(Crystal *cryst, double m); + +#endif /* CRYSTAL_H */ diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c index 160df719..e6ff36b1 100644 --- a/libcrystfel/src/dirax.c +++ b/libcrystfel/src/dirax.c @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -53,6 +53,7 @@ #include "dirax.h" #include "utils.h" #include "peaks.h" +#include "cell-utils.h" #define DIRAX_VERBOSE 0 @@ -68,6 +69,13 @@ typedef enum { } DirAxInputType; +struct dirax_private { + IndexingMethod indm; + float *ltl; + UnitCell *template; +}; + + struct dirax_data { /* DirAx auto-indexing low-level stuff */ @@ -83,12 +91,59 @@ struct dirax_data { int read_cell; int best_acl; int best_acl_nh; - int acls_tried[MAX_CELL_CANDIDATES]; + int acls_tried[MAX_DIRAX_CELL_CANDIDATES]; int n_acls_tried; + UnitCell *cur_cell; + int done; + int success; + + struct dirax_private *dp; }; +static int check_cell(struct dirax_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + if ( dp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 1); + if ( out == NULL ) return 0; + + } else if ( dp->indm & INDEXING_CHECK_CELL_AXES ) { + + out = match_cell(cell, dp->template, 0, dp->ltl, 0); + if ( out == NULL ) return 0; + + } else { + out = cell_new_from_cell(cell); + } + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + cell_free(out); + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + + static void dirax_parseline(const char *line, struct image *image, struct dirax_data *dirax) { @@ -119,7 +174,7 @@ static void dirax_parseline(const char *line, struct image *image, if ( line[i] == 'R' ) rf = 1; if ( (line[i] == 'D') && rf ) { dirax->read_cell = 1; - image->candidate_cells[image->ncells] = cell_new(); + dirax->cur_cell = cell_new(); return; } i++; @@ -127,6 +182,7 @@ static void dirax_parseline(const char *line, struct image *image, /* Parse unit cell vectors as appropriate */ if ( dirax->read_cell == 1 ) { + /* First row of unit cell values */ float ax, ay, az; int r; @@ -136,14 +192,16 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_a(image->candidate_cells[image->ncells], + cell_set_cartesian_a(dirax->cur_cell, ax*1e-10, ay*1e-10, az*1e-10); dirax->read_cell++; return; + } else if ( dirax->read_cell == 2 ) { + /* Second row of unit cell values */ float bx, by, bz; int r; @@ -153,14 +211,16 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_b(image->candidate_cells[image->ncells], + cell_set_cartesian_b(dirax->cur_cell, bx*1e-10, by*1e-10, bz*1e-10); dirax->read_cell++; return; + } else if ( dirax->read_cell == 3 ) { + /* Third row of unit cell values */ float cx, cy, cz; int r; @@ -170,13 +230,21 @@ static void dirax_parseline(const char *line, struct image *image, ERROR("Couldn't understand cell line:\n"); ERROR("'%s'\n", line); dirax->read_cell = 0; - cell_free(image->candidate_cells[image->ncells]); + cell_free(dirax->cur_cell); return; } - cell_set_cartesian_c(image->candidate_cells[image->ncells++], + cell_set_cartesian_c(dirax->cur_cell, cx*1e-10, cy*1e-10, cz*1e-10); dirax->read_cell = 0; + + /* Finished reading a cell. Time to check it... */ + if ( check_cell(dirax->dp, image, dirax->cur_cell) ) { + dirax->done = 1; + dirax->success = 1; + } + return; + } dirax->read_cell = 0; @@ -351,8 +419,7 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax) switch ( type ) { - case DIRAX_INPUT_LINE : - + case DIRAX_INPUT_LINE : block_buffer = malloc(i+1); memcpy(block_buffer, dirax->rbuffer, i); block_buffer[i] = '\0'; @@ -366,20 +433,17 @@ static int dirax_readable(struct image *image, struct dirax_data *dirax) endbit_length = i+2; break; - case DIRAX_INPUT_PROMPT : - + case DIRAX_INPUT_PROMPT : dirax_send_next(image, dirax); endbit_length = i+7; break; - case DIRAX_INPUT_ACL : - + case DIRAX_INPUT_ACL : dirax_send_next(image,dirax ); endbit_length = i+10; break; - default : - + default : /* Obviously, this never happens :) */ ERROR("Unrecognised DirAx input mode! " "I don't know how to understand DirAx\n"); @@ -456,7 +520,7 @@ static void write_drx(struct image *image) } -void run_dirax(struct image *image) +int run_dirax(struct image *image, IndexingPrivate *ipriv) { unsigned int opts; int status; @@ -468,13 +532,13 @@ void run_dirax(struct image *image) dirax = malloc(sizeof(struct dirax_data)); if ( dirax == NULL ) { ERROR("Couldn't allocate memory for DirAx data.\n"); - return; + return 0; } dirax->pid = forkpty(&dirax->pty, NULL, NULL, NULL); if ( dirax->pid == -1 ) { ERROR("Failed to fork for DirAx: %s\n", strerror(errno)); - return; + return 0; } if ( dirax->pid == 0 ) { @@ -505,6 +569,9 @@ void run_dirax(struct image *image) dirax->read_cell = 0; dirax->n_acls_tried = 0; dirax->best_acl_nh = 0; + dirax->done = 0; + dirax->success = 0; + dirax->dp = (struct dirax_private *)ipriv; do { @@ -543,7 +610,7 @@ void run_dirax(struct image *image) rval = 1; } - } while ( !rval ); + } while ( !rval && !dirax->success ); close(dirax->pty); free(dirax->rbuffer); @@ -553,5 +620,45 @@ void run_dirax(struct image *image) ERROR("DirAx doesn't seem to be working properly.\n"); } + rval = dirax->success; free(dirax); + return rval; +} + + +IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) +{ + struct dirax_private *dp; + int need_cell = 0; + + if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; + if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; + + if ( need_cell && (cell == NULL) ) { + ERROR("DirAx needs a unit cell for this set of flags.\n"); + return NULL; + } + + /* Flags that DirAx knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS + | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS; + + dp = malloc(sizeof(struct dirax_private)); + if ( dp == NULL ) return NULL; + + dp->ltl = ltl; + dp->template = cell; + dp->indm = *indm; + + return (IndexingPrivate *)dp; +} + + +void dirax_cleanup(IndexingPrivate *pp) +{ + struct dirax_private *p; + p = (struct dirax_private *)pp; + free(p); } diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h index ce3cd4b0..6be8451a 100644 --- a/libcrystfel/src/dirax.h +++ b/libcrystfel/src/dirax.h @@ -3,11 +3,11 @@ * * Invoke the DirAx auto-indexing program * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010,2012 Thomas White <taw@physics.org> + * 2010,2012-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -33,10 +33,15 @@ #include <config.h> #endif -#include "utils.h" +#include "index.h" +extern int run_dirax(struct image *image, IndexingPrivate *ipriv); -extern void run_dirax(struct image *image); +extern IndexingPrivate *dirax_prepare(IndexingMethod *indm, + UnitCell *cell, const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); +extern void dirax_cleanup(IndexingPrivate *pp); #endif /* DIRAX_H */ diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 60fab488..704baa51 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -3,11 +3,11 @@ * * Geometry of diffraction * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -33,6 +33,7 @@ #include <stdlib.h> #include <assert.h> +#include <fenv.h> #include "utils.h" #include "cell.h" @@ -43,6 +44,7 @@ #include "reflist.h" #include "reflist-utils.h" #include "symmetry.h" +#include "geometry.h" static signed int locate_peak(double x, double y, double z, double k, @@ -117,7 +119,7 @@ static double partiality(double rlow, double rhigh, double r) } -static Reflection *check_reflection(struct image *image, +static Reflection *check_reflection(struct image *image, Crystal *cryst, signed int h, signed int k, signed int l, double xl, double yl, double zl) { @@ -131,6 +133,9 @@ static Reflection *check_reflection(struct image *image, double klow, khigh; /* Wavenumber */ Reflection *refl; double cet, cez; + double pr; + + pr = crystal_get_profile_radius(cryst); /* "low" gives the largest Ewald sphere (wavelength short => k large) * "high" gives the smallest Ewald sphere (wavelength long => k small) @@ -152,8 +157,8 @@ static Reflection *check_reflection(struct image *image, rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ if ( (signbit(rlow) == signbit(rhigh)) - && (fabs(rlow) > image->profile_radius) - && (fabs(rhigh) > image->profile_radius) ) return NULL; + && (fabs(rlow) > pr) + && (fabs(rhigh) > pr) ) return NULL; /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the @@ -164,26 +169,26 @@ static Reflection *check_reflection(struct image *image, * et al. (1979). */ clamp_low = 0; clamp_high = 0; - if ( rlow < -image->profile_radius ) { - rlow = -image->profile_radius; + if ( rlow < -pr ) { + rlow = -pr; clamp_low = -1; } - if ( rlow > +image->profile_radius ) { - rlow = +image->profile_radius; + if ( rlow > +pr ) { + rlow = +pr; clamp_low = +1; } - if ( rhigh < -image->profile_radius ) { - rhigh = -image->profile_radius; + if ( rhigh < -pr ) { + rhigh = -pr; clamp_high = -1; } - if ( rhigh > +image->profile_radius ) { - rhigh = +image->profile_radius; + if ( rhigh > +pr ) { + rhigh = +pr; clamp_high = +1; } assert(clamp_low >= clamp_high); /* Calculate partiality */ - part = partiality(rlow, rhigh, image->profile_radius); + part = partiality(rlow, rhigh, pr); /* Locate peak on detector. */ p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det, &xda, &yda); @@ -205,7 +210,7 @@ static Reflection *check_reflection(struct image *image, } -RefList *find_intersections(struct image *image, UnitCell *cell) +RefList *find_intersections(struct image *image, Crystal *cryst) { double ax, ay, az; double bx, by, bz; @@ -217,6 +222,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell) int hmax, kmax, lmax; double mres; signed int h, k, l; + UnitCell *cell; + + cell = crystal_get_cell(cryst); + if ( cell == NULL ) return NULL; reflections = reflist_new(); @@ -262,7 +271,7 @@ RefList *find_intersections(struct image *image, UnitCell *cell) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - refl = check_reflection(image, h, k, l, xl, yl, zl); + refl = check_reflection(image, cryst, h, k, l, xl, yl, zl); if ( refl != NULL ) { add_refl_to_list(refl, reflections); @@ -276,8 +285,68 @@ RefList *find_intersections(struct image *image, UnitCell *cell) } +RefList *select_intersections(struct image *image, Crystal *cryst) +{ + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + const double min_dist = 0.25; + RefList *list; + int i; + + /* Round towards nearest */ + fesetround(1); + + /* Cell basis vectors for this image */ + cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az, + &bx, &by, &bz, &cx, &cy, &cz); + + list = reflist_new(); + if ( list == NULL ) return NULL; + + /* Loop over peaks, checking proximity to nearest reflection */ + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f; + struct rvec q; + double h, k, l, hd, kd, ld; + double dsq; + + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + /* Reciprocal space position of found peak */ + q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); + + /* Decimal and fractional Miller indices of nearest + * reciprocal lattice point */ + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + /* Check distance */ + dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); + + if ( sqrt(dsq) < min_dist ) { + + Reflection *refl; + + refl = add_refl(list, h, k, l); + set_detector_pos(refl, sqrt(dsq), f->fs, f->ss); + + } + + } + + return list; +} + + /* Calculate partialities and apply them to the image's reflections */ -void update_partialities(struct image *image) +void update_partialities(Crystal *cryst) { Reflection *refl; RefListIterator *iter; @@ -285,14 +354,15 @@ void update_partialities(struct image *image) double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; + struct image *image = crystal_get_image(cryst); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); /* Scratch list to give check_reflection() something to add to */ predicted = reflist_new(); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -309,7 +379,7 @@ void update_partialities(struct image *image) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - vals = check_reflection(image, h, k, l, xl, yl, zl); + vals = check_reflection(image, cryst, h, k, l, xl, yl, zl); if ( vals == NULL ) { set_redundancy(refl, 0); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 94b496e6..aecdc28a 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -3,12 +3,12 @@ * * Geometry of diffraction * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2012 Richard Kirian * * This file is part of CrystFEL. @@ -43,9 +43,10 @@ extern "C" { #endif -RefList *find_intersections(struct image *image, UnitCell *cell); +extern RefList *find_intersections(struct image *image, Crystal *cryst); +extern RefList *select_intersections(struct image *image, Crystal *cryst); -void update_partialities(struct image *image); +extern void update_partialities(Crystal *cryst); #ifdef __cplusplus } diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 093e20f2..6e4fd40a 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -159,3 +159,21 @@ void image_remove_feature(ImageFeatureList *flist, int idx) { flist->features[idx].valid = 0; } + + +void image_add_crystal(struct image *image, Crystal *cryst) +{ + Crystal **crs; + int n; + + n = image->n_crystals; + crs = realloc(image->crystals, (n+1)*sizeof(Crystal *)); + if ( crs == NULL ) { + ERROR("Failed to allocate memory for crystals.\n"); + return; + } + + crs[n] = cryst; + image->crystals = crs; + image->n_crystals = n+1; +} diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index afa9e4a7..3739c01f 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -3,11 +3,11 @@ * * Handle images and image features * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2009-2012 Thomas White <taw@physics.org> + * 2009-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -41,9 +41,8 @@ #include "cell.h" #include "detector.h" #include "reflist.h" - - -#define MAX_CELL_CANDIDATES (32) +#include "crystal.h" +#include "index.h" /* Structure describing a feature in an image */ @@ -79,10 +78,10 @@ typedef struct _imagefeaturelist ImageFeatureList; * uint16_t *flags; * double *twotheta; * - * UnitCell *indexed_cell; - * UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; - * int ncells; - + * Crystal **crystals; + * int n_crystals; + * IndexingMethod indexed_by; + * * struct detector *det; * struct beam_params *beam; * char *filename; @@ -90,22 +89,15 @@ typedef struct _imagefeaturelist ImageFeatureList; * * int id; * - * double m; - * * double lambda; * double div; * double bw; - * double osf; - * double profile_radius; - * int pr_dud; - * double diffracting_resolution; * * int width; * int height; * - * RefList *reflections; - * long long unsigned int n_saturated; - * + * long long int num_peaks; + * long long int num_saturated_peaks; * ImageFeatureList *features; * }; * </programlisting> @@ -121,14 +113,9 @@ typedef struct _imagefeaturelist ImageFeatureList; * by-product of the scattering vector calculation and can be used later for * calculating intensities from differential scattering cross sections. * - * <structfield>candidate_cells</structfield> is an array of unit cells directly - * returned by the low-level indexing system. <structfield>ncells</structfield> - * is the number of candidate unit cells which were found. The maximum number - * of cells which may be returned is <function>MAX_CELL_CANDIDATES</function>. - * <structfield>indexed_cell</structfield> contains the "correct" unit cell - * after cell reduction or matching has been performed. The job of the cell - * reduction is to convert the list of candidate cells into a single indexed - * cell, or <function>NULL</function> on failure. + * <structfield>crystals</structfield> is an array of %Crystal directly + * returned by the low-level indexing system. <structfield>n_crystals</structfield> + * is the number of crystals which were found in the image. * * <structfield>copyme</structfield> represents a list of HDF5 fields to copy * to the output stream. @@ -141,9 +128,9 @@ struct image { uint16_t *flags; double *twotheta; - UnitCell *indexed_cell; - UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; - int ncells; + Crystal **crystals; + int n_crystals; + IndexingMethod indexed_by; struct detector *det; struct beam_params *beam; /* The nominal beam parameters */ @@ -153,26 +140,17 @@ struct image { int id; /* ID number of the thread * handling this image */ - /* Information about the crystal */ - double m; /* Mosaicity in radians */ - /* Per-shot radiation values */ double lambda; /* Wavelength in m */ double div; /* Divergence in radians */ double bw; /* Bandwidth as a fraction */ - double osf; /* Overall scaling factor */ - double profile_radius; /* Radius of reflection */ - int pr_dud; /* Post refinement failed */ - double diffracting_resolution; /* Max 1/d in m^-1 */ int width; int height; - /* Integrated (or about-to-be-integrated) reflections */ - RefList *reflections; - long long int n_saturated; /* Number of overloads */ - /* Detected peaks */ + long long int num_peaks; + long long int num_saturated_peaks; ImageFeatureList *features; }; @@ -196,4 +174,6 @@ extern struct imagefeature *image_feature_closest(ImageFeatureList *flist, extern int image_feature_count(ImageFeatureList *flist); extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx); +extern void image_add_crystal(struct image *image, Crystal *cryst); + #endif /* IMAGE_H */ diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 86ab8c1a..89fdec4f 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -3,12 +3,12 @@ * * Perform indexing (somehow) * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2010-2011 Richard Kirian <rkirian@asu.edu> * 2012 Lorenzo Galli * @@ -53,49 +53,34 @@ #include "cell-utils.h" -/* Base class constructor for unspecialised indexing private data */ -IndexingPrivate *indexing_private(IndexingMethod indm) -{ - struct _indexingprivate *priv; - priv = calloc(1, sizeof(struct _indexingprivate)); - priv->indm = indm; - return priv; -} - - -static const char *maybes(int n) -{ - if ( n == 1 ) return ""; - return "s"; -} - - IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, const char *filename, struct detector *det, - double nominal_photon_energy) + struct beam_params *beam, float *ltl) { int n; int nm = 0; IndexingPrivate **iprivs; while ( indm[nm] != INDEXING_NONE ) nm++; - STATUS("Preparing %i indexing method%s.\n", nm, maybes(nm)); iprivs = malloc((nm+1) * sizeof(IndexingPrivate *)); for ( n=0; n<nm; n++ ) { - switch ( indm[n] ) { + int i; + IndexingMethod in; - case INDEXING_NONE : - ERROR("Tried to prepare INDEXING_NONE!\n"); - break; + in = indm[n]; + + switch ( indm[n] & INDEXING_METHOD_MASK ) { case INDEXING_DIRAX : - iprivs[n] = indexing_private(indm[n]); + iprivs[n] = dirax_prepare(&indm[n], cell, filename, + det, beam, ltl); break; case INDEXING_MOSFLM : - iprivs[n] = indexing_private(indm[n]); + iprivs[n] = mosflm_prepare(&indm[n], cell, filename, + det, beam, ltl); break; case INDEXING_GRAINSPOTTER : @@ -103,11 +88,38 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, break; case INDEXING_REAX : - iprivs[n] = reax_prepare(); + iprivs[n] = reax_prepare(&indm[n], cell, filename, + det, beam, ltl); + break; + + default : + ERROR("Don't know how to prepare indexing method %i\n", + indm[n]); break; } + if ( iprivs[n] == NULL ) return NULL; + + STATUS("Prepared indexing method %i: %s\n", + n, indexer_str(indm[n])); + + if ( in != indm[n] ) { + ERROR("Note: flags were altered to take into account " + "the limitations of the indexing method.\n"); + } + + for ( i=0; i<n; i++ ) { + if ( indm[i] == indm[n] ) { + ERROR("Duplicate indexing method.\n"); + ERROR("Have you specified some flags which " + "aren't accepted by one of your " + "chosen indexing methods?\n"); + return NULL; + } + } + + } iprivs[n] = NULL; @@ -115,26 +127,26 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, } -void cleanup_indexing(IndexingPrivate **priv) +void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs) { int n = 0; - if ( priv == NULL ) return; /* Nothing to do */ + if ( indms == NULL ) return; /* Nothing to do */ + if ( privs == NULL ) return; /* Nothing to do */ - while ( priv[n] != NULL ) { + while ( indms[n] != INDEXING_NONE ) { - switch ( priv[n]->indm ) { + switch ( indms[n] & INDEXING_METHOD_MASK ) { case INDEXING_NONE : - free(priv[n]); break; case INDEXING_DIRAX : - free(priv[n]); + dirax_cleanup(privs[n]); break; case INDEXING_MOSFLM : - free(priv[n]); + mosflm_cleanup(privs[n]); break; case INDEXING_GRAINSPOTTER : @@ -142,7 +154,12 @@ void cleanup_indexing(IndexingPrivate **priv) break; case INDEXING_REAX : - reax_cleanup(priv[n]); + reax_cleanup(privs[n]); + break; + + default : + ERROR("Don't know how to clean up indexing method %i\n", + indms[n]); break; } @@ -175,146 +192,213 @@ void map_all_peaks(struct image *image) } -void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, - int cellr, int verbose, IndexingPrivate **ipriv, - int config_insane, const float *ltl) +/* Return non-zero for "success" */ +static int try_indexer(struct image *image, IndexingMethod indm, + IndexingPrivate *ipriv) +{ + switch ( indm & INDEXING_METHOD_MASK ) { + + case INDEXING_NONE : + break; + + case INDEXING_DIRAX : + return run_dirax(image, ipriv); + break; + + case INDEXING_MOSFLM : + return run_mosflm(image, ipriv); + break; + + case INDEXING_REAX : + return reax_index(image, ipriv); + break; + + default : + ERROR("Unrecognised indexing method: %i\n", indm); + break; + + } + + return 0; +} + + +void index_pattern(struct image *image, + IndexingMethod *indms, IndexingPrivate **iprivs) { - int i; int n = 0; - if ( indm == NULL ) return; + if ( indms == NULL ) return; map_all_peaks(image); - image->indexed_cell = NULL; + image->crystals = NULL; + image->n_crystals = 0; - while ( indm[n] != INDEXING_NONE ) { + while ( indms[n] != INDEXING_NONE ) { - image->ncells = 0; + if ( try_indexer(image, indms[n], iprivs[n]) ) break; + n++; - /* Index as appropriate */ - switch ( indm[n] ) { + } - case INDEXING_NONE : - return; + image->indexed_by = indms[n]; +} - case INDEXING_DIRAX : - run_dirax(image); - break; - case INDEXING_MOSFLM : - run_mosflm(image, cell); - break; +/* Set the indexer flags for "raw mode" ("--cell-reduction=none") */ +static IndexingMethod set_raw(IndexingMethod a) +{ + /* Disable all unit cell checks */ + a &= ~(INDEXING_CHECK_CELL_COMBINATIONS | INDEXING_CHECK_CELL_AXES); + return a; +} - case INDEXING_GRAINSPOTTER : - run_grainspotter(image, cell); - break; - case INDEXING_REAX : - reax_index(ipriv[n], image, cell); - break; +/* Set the indexer flags for "bad mode" ("--insane) */ +static IndexingMethod set_bad(IndexingMethod a) +{ + /* Disable the peak check */ + return a & ~INDEXING_CHECK_PEAKS; +} - } - if ( image->ncells == 0 ) { - n++; - continue; - } - for ( i=0; i<image->ncells; i++ ) { +/* Set the indexer flags for "axes mode" ("--cell-reduction=compare") */ +static IndexingMethod set_axes(IndexingMethod a) +{ + return (a & ~INDEXING_CHECK_CELL_COMBINATIONS) + | INDEXING_CHECK_CELL_AXES; +} - UnitCell *new_cell = NULL; - UnitCell *cand = image->candidate_cells[i]; - if ( verbose ) { - STATUS("--------------------\n"); - STATUS("Candidate cell %i (before matching):\n", - i); - cell_print(image->candidate_cells[i]); - STATUS("--------------------\n"); - } +/* Set the indexer flags for "combination mode" ("--cell-reduction=reduce") */ +static IndexingMethod set_comb(IndexingMethod a) +{ + return (a & ~INDEXING_CHECK_CELL_AXES) + | INDEXING_CHECK_CELL_COMBINATIONS; +} - /* Match or reduce the cell as appropriate */ - switch ( cellr ) { - case CELLR_NONE : - new_cell = cell_new_from_cell(cand); - break; +/* Set the indexer flags for "use no lattice type information" */ +static IndexingMethod set_nolattice(IndexingMethod a) +{ + return a & ~INDEXING_USE_LATTICE_TYPE; +} - case CELLR_REDUCE : - new_cell = match_cell(cand, cell, verbose, - ltl, 1); - break; - case CELLR_COMPARE : - new_cell = match_cell(cand, cell, verbose, - ltl, 0); - break; +/* Set the indexer flags for "use lattice type information" */ +static IndexingMethod set_lattice(IndexingMethod a) +{ + return a | INDEXING_USE_LATTICE_TYPE; +} - case CELLR_COMPARE_AB : - new_cell = match_cell_ab(cand, cell); - break; - } +char *indexer_str(IndexingMethod indm) +{ + char *str; - /* No cell? Move on to the next candidate */ - if ( new_cell == NULL ) continue; + str = malloc(32); + if ( str == NULL ) { + ERROR("Failed to allocate string.\n"); + return NULL; + } + str[0] = '\0'; - /* Sanity check */ - image->indexed_cell = new_cell; - if ( !config_insane && !peak_sanity_check(image) ) { - cell_free(new_cell); - image->indexed_cell = NULL; - continue; - } + switch ( indm & INDEXING_METHOD_MASK ) { - goto done; /* Success */ + case INDEXING_NONE : + strcpy(str, "none"); + return str; - } + case INDEXING_DIRAX : + strcpy(str, "dirax"); + break; - for ( i=0; i<image->ncells; i++ ) { - cell_free(image->candidate_cells[i]); - image->candidate_cells[i] = NULL; - } + case INDEXING_MOSFLM : + strcpy(str, "mosflm"); + break; - /* Move on to the next indexing method */ - n++; + case INDEXING_REAX : + strcpy(str, "reax"); + break; + + default : + ERROR("Unrecognised indexing method %i\n", indm); + strcpy(str, "(unknown)"); + break; } -done: - for ( i=0; i<image->ncells; i++ ) { - /* May free(NULL) if all algorithms were tried and no success */ - cell_free(image->candidate_cells[i]); + if ( indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + strcat(str, "-comb"); + } else if ( indm & INDEXING_CHECK_CELL_AXES ) { + strcat(str, "-axes"); + } else { + strcat(str, "-raw"); } + + if ( !(indm & INDEXING_CHECK_PEAKS) ) { + strcat(str, "-bad"); + } + + if ( indm & INDEXING_USE_LATTICE_TYPE ) { + strcat(str, "-latt"); + } else { + strcat(str, "-nolatt"); + } + + return str; } -IndexingMethod *build_indexer_list(const char *str, int *need_cell) +IndexingMethod *build_indexer_list(const char *str) { int n, i; char **methods; IndexingMethod *list; - int tmp; - - if ( need_cell == NULL ) need_cell = &tmp; - *need_cell = 0; + int nmeth = 0; - n = assplode(str, ",", &methods, ASSPLODE_NONE); + n = assplode(str, ",-", &methods, ASSPLODE_NONE); list = malloc((n+1)*sizeof(IndexingMethod)); + nmeth = -1; /* So that the first method is #0 */ for ( i=0; i<n; i++ ) { if ( strcmp(methods[i], "dirax") == 0) { - list[i] = INDEXING_DIRAX; + list[++nmeth] = INDEXING_DEFAULTS_DIRAX; + } else if ( strcmp(methods[i], "mosflm") == 0) { - list[i] = INDEXING_MOSFLM; + list[++nmeth] = INDEXING_DEFAULTS_MOSFLM; + } else if ( strcmp(methods[i], "grainspotter") == 0) { - list[i] = INDEXING_GRAINSPOTTER; + list[++nmeth] = INDEXING_DEFAULTS_GRAINSPOTTER; + } else if ( strcmp(methods[i], "reax") == 0) { - list[i] = INDEXING_REAX; - *need_cell = 1; + list[++nmeth] = INDEXING_DEFAULTS_REAX; + + } else if ( strcmp(methods[i], "none") == 0) { + list[++nmeth] = INDEXING_NONE; + return list; + + } else if ( strcmp(methods[i], "raw") == 0) { + list[nmeth] = set_raw(list[nmeth]); + + } else if ( strcmp(methods[i], "bad") == 0) { + list[nmeth] = set_bad(list[nmeth]); + + } else if ( strcmp(methods[i], "comb") == 0) { + list[nmeth] = set_comb(list[nmeth]); /* Default */ + + } else if ( strcmp(methods[i], "axes") == 0) { + list[nmeth] = set_axes(list[nmeth]); + + } else if ( strcmp(methods[i], "latt") == 0) { + list[nmeth] = set_lattice(list[nmeth]); + + } else if ( strcmp(methods[i], "nolatt") == 0) { + list[nmeth] = set_nolattice(list[nmeth]); + } else { - ERROR("Unrecognised indexing method '%s'\n", - methods[i]); + ERROR("Bad list of indexing methods: '%s'\n", str); return NULL; } @@ -322,7 +406,7 @@ IndexingMethod *build_indexer_list(const char *str, int *need_cell) } free(methods); - list[i] = INDEXING_NONE; + list[++nmeth] = INDEXING_NONE; return list; } diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 2df7a311..a47ee9a0 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -3,13 +3,13 @@ * * Perform indexing (somehow) * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2010 Richard Kirian * 2012 Lorenzo Galli * @@ -38,11 +38,22 @@ #endif +#include "beam-parameters.h" #include "cell.h" #include "image.h" #include "detector.h" +#define INDEXING_DEFAULTS_DIRAX (INDEXING_DIRAX | INDEXING_CHECK_PEAKS \ + | INDEXING_CHECK_CELL_COMBINATIONS) + +#define INDEXING_DEFAULTS_MOSFLM (INDEXING_MOSFLM | INDEXING_CHECK_PEAKS \ + | INDEXING_CHECK_CELL_COMBINATIONS \ + | INDEXING_USE_LATTICE_TYPE) + +#define INDEXING_DEFAULTS_REAX (INDEXING_REAX | INDEXING_USE_LATTICE_TYPE \ + | INDEXING_CHECK_PEAKS) + /** * IndexingMethod: * @INDEXING_NONE: No indexing to be performed @@ -53,21 +64,28 @@ * An enumeration of all the available indexing methods. **/ typedef enum { - INDEXING_NONE, - INDEXING_DIRAX, - INDEXING_MOSFLM, - INDEXING_GRAINSPOTTER, - INDEXING_REAX, -} IndexingMethod; -/* Cell reduction methods */ -enum { - CELLR_NONE, - CELLR_REDUCE, - CELLR_COMPARE, - CELLR_COMPARE_AB, -}; + INDEXING_NONE = 0, + + /* The core indexing methods themselves */ + INDEXING_DIRAX = 1, + INDEXING_MOSFLM = 2, + INDEXING_REAX = 3, + INDEXING_GRAINSPOTTER = 4, + + /* Bits at the top of the IndexingMethod are flags which modify the + * behaviour of the indexer, at the moment just by adding checks. */ + INDEXING_CHECK_CELL_COMBINATIONS = 256, + INDEXING_CHECK_CELL_AXES = 512, + INDEXING_CHECK_PEAKS = 1024, + INDEXING_USE_LATTICE_TYPE = 2048, + +} IndexingMethod; + +/* This defines the bits in "IndexingMethod" which are used to represent the + * core of the indexing method */ +#define INDEXING_METHOD_MASK (0xff) /** @@ -76,24 +94,19 @@ enum { * This is an opaque data structure containing information needed by the * indexing method. **/ -typedef struct _indexingprivate IndexingPrivate; +typedef void *IndexingPrivate; -extern IndexingPrivate *indexing_private(IndexingMethod indm); +extern IndexingMethod *build_indexer_list(const char *str); +extern char *indexer_str(IndexingMethod indm); extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell, - const char *filename, - struct detector *det, - double nominal_photon_energy); - -extern void map_all_peaks(struct image *image); - -extern void index_pattern(struct image *image, UnitCell *cell, - IndexingMethod *indm, int cellr, int verbose, - IndexingPrivate **priv, int config_insane, - const float *ltl); + const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); -extern void cleanup_indexing(IndexingPrivate **priv); +extern void index_pattern(struct image *image, + IndexingMethod *indms, IndexingPrivate **iprivs); -extern IndexingMethod *build_indexer_list(const char *str, int *need_cell); +extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs); #endif /* INDEX_H */ diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 548e3834..a03e621d 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian <rkirian@asu.edu> - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -80,6 +80,7 @@ #include "mosflm.h" #include "utils.h" #include "peaks.h" +#include "cell-utils.h" #define MOSFLM_VERBOSE 0 @@ -92,6 +93,14 @@ typedef enum { } MOSFLMInputType; + +struct mosflm_private { + IndexingMethod indm; + float *ltl; + UnitCell *template; +}; + + struct mosflm_data { /* MOSFLM auto-indexing low-level stuff */ @@ -107,10 +116,85 @@ struct mosflm_data { char sptfile[128]; int step; int finished_ok; - UnitCell *target_cell; + int done; + int success; + + struct mosflm_private *mp; }; +static int check_cell(struct mosflm_private *mp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + /* If we sent lattice information, make sure that we got back what we + * asked for, not (e.g.) some "H" version of a rhombohedral R cell */ + if ( mp->indm & INDEXING_USE_LATTICE_TYPE ) { + + LatticeType latt_m, latt_r; + char cen_m, cen_r; + + /* What we asked for */ + latt_r = cell_get_lattice_type(mp->template); + cen_r = cell_get_centering(mp->template); + + /* What we got back */ + latt_m = cell_get_lattice_type(cell); + cen_m = cell_get_centering(cell); + + if ( latt_r != latt_m ) { + ERROR("Lattice type produced by MOSFLM (%i) does not " + "match what was requested (%i). " + "Please report this.\n", latt_m, latt_r); + return 0; + } + + if ( cen_r != cen_m ) { + ERROR("Centering produced by MOSFLM (%c) does not " + "match what was requested (%c). " + "Please report this.\n", cen_m, cen_r); + return 0; + } + + } + + if ( mp->indm & INDEXING_CHECK_CELL_COMBINATIONS ) { + + out = match_cell(cell, mp->template, 0, mp->ltl, 1); + if ( out == NULL ) return 0; + + } else if ( mp->indm & INDEXING_CHECK_CELL_AXES ) { + + out = match_cell(cell, mp->template, 0, mp->ltl, 0); + if ( out == NULL ) return 0; + + } else { + out = cell_new_from_cell(cell); + } + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + if ( mp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + cell_free(out); + crystal_free(cr); + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + static void mosflm_parseline(const char *line, struct image *image, struct mosflm_data *dirax) @@ -130,14 +214,54 @@ static void mosflm_parseline(const char *line, struct image *image, } -static int read_newmat(const char *filename, struct image *image) +/* This is the opposite of spacegroup_for_lattice() below. */ +static LatticeType spacegroup_to_lattice(const char *sg) +{ + LatticeType latt; + + if ( sg[1] == '1' ) { + latt = L_TRICLINIC; + } else if ( strncmp(sg+1, "23", 2) == 0 ) { + latt = L_CUBIC; + } else if ( strncmp(sg+1, "222", 3) == 0 ) { + latt = L_ORTHORHOMBIC; + } else if ( sg[1] == '2' ) { + latt = L_MONOCLINIC; + } else if ( sg[1] == '4' ) { + latt = L_TETRAGONAL; + } else if ( sg[1] == '6' ) { + latt = L_HEXAGONAL; + } else if ( sg[1] == '3' ) { + if ( sg[0] == 'H' ) { + latt = L_HEXAGONAL; + } else { + latt = L_RHOMBOHEDRAL; + } + } else { + ERROR("Couldn't understand '%s'\n", sg); + latt = L_TRICLINIC; + } + + return latt; +} + + + +static int read_newmat(struct mosflm_data *mosflm, const char *filename, + struct image *image) { - FILE * fh; + FILE *fh; float asx, asy, asz; float bsx, bsy, bsz; float csx, csy, csz; int n; double c; + UnitCell *cell; + char symm[32]; + char *rval; + int i; + char cen; + LatticeType latt; fh = fopen(filename, "r"); if ( fh == NULL ) { @@ -150,23 +274,54 @@ static int read_newmat(const char *filename, struct image *image) STATUS("Fewer than 9 parameters found in NEWMAT file.\n"); return 1; } + + /* Skip the next six lines */ + for ( i=0; i<6; i++ ) { + char tmp[1024]; + rval = fgets(tmp, 1024, fh); + if ( rval == NULL ) { + ERROR("Failed to read newmat file.\n"); + return 1; + } + } + + rval = fgets(symm, 32, fh); + if ( rval == NULL ) { + ERROR("Failed to read newmat file.\n"); + return 1; + } + fclose(fh); + if ( strncmp(symm, "SYMM ", 5) != 0 ) { + ERROR("Bad 'SYMM' line from MOSFLM.\n"); + return 1; + } + cen = symm[5]; + latt = spacegroup_to_lattice(symm+5); + /* MOSFLM "A" matrix is multiplied by lambda, so fix this */ c = 1.0/image->lambda; - image->candidate_cells[0] = cell_new(); + cell = cell_new(); /* The relationship between the coordinates in the spot file and the * resulting matrix is diabolically complicated. This transformation * seems to work, but is not derived by working through all the * transformations. */ - cell_set_reciprocal(image->candidate_cells[0], + cell_set_reciprocal(cell, -asy*c, -asz*c, asx*c, -bsy*c, -bsz*c, bsx*c, -csy*c, -csz*c, csx*c); + cell_set_centering(cell, cen); + cell_set_lattice_type(cell, latt); + + if ( check_cell(mosflm->mp, image, cell) ) { + mosflm->success = 1; + mosflm->done = 1; + } - image->ncells = 1; + cell_free(cell); return 0; } @@ -319,7 +474,7 @@ static const char *spacegroup_for_lattice(UnitCell *cell) if ( centering != 'H' ) { g = "6"; } else { - g = "32"; + g = "3"; } break; @@ -343,8 +498,8 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) char tmp[256]; double wavelength; - switch ( mosflm->step ) { - + switch ( mosflm->step ) + { case 1 : mosflm_sendline("DETECTOR ROTATION HORIZONTAL" " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL" @@ -352,13 +507,22 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) break; case 2 : - if ( mosflm->target_cell != NULL ) { + if ( (mosflm->mp->indm & INDEXING_USE_LATTICE_TYPE) + && (mosflm->mp->template != NULL) ) + { const char *symm; - symm = spacegroup_for_lattice(mosflm->target_cell); + + if ( cell_get_lattice_type(mosflm->mp->template) + == L_RHOMBOHEDRAL ) { + mosflm_sendline("CRYSTAL R\n", mosflm); + } + + symm = spacegroup_for_lattice(mosflm->mp->template); snprintf(tmp, 255, "SYMM %s\n", symm); mosflm_sendline(tmp, mosflm); + } else { - mosflm_sendline("SYMM P1\n", mosflm); + mosflm_sendline("\n", mosflm); } break; @@ -403,10 +567,9 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm) mosflm->finished_ok = 1; break; - default: + default: mosflm_sendline("exit\n", mosflm); return; - } mosflm->step++; @@ -460,8 +623,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) switch ( type ) { - case MOSFLM_INPUT_LINE : - + case MOSFLM_INPUT_LINE : block_buffer = malloc(i+1); memcpy(block_buffer, mosflm->rbuffer, i); block_buffer[i] = '\0'; @@ -475,12 +637,12 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) endbit_length = i+2; break; - case MOSFLM_INPUT_PROMPT : + case MOSFLM_INPUT_PROMPT : mosflm_send_next(image, mosflm); endbit_length = i+7; break; - default : + default : /* Obviously, this never happens :) */ ERROR("Unrecognised MOSFLM input mode! " @@ -527,7 +689,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) } -void run_mosflm(struct image *image, UnitCell *cell) +int run_mosflm(struct image *image, IndexingPrivate *ipriv) { struct mosflm_data *mosflm; unsigned int opts; @@ -537,11 +699,9 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm = malloc(sizeof(struct mosflm_data)); if ( mosflm == NULL ) { ERROR("Couldn't allocate memory for MOSFLM data.\n"); - return; + return 0; } - mosflm->target_cell = cell; - snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id); write_img(image, mosflm->imagefile); /* Dummy image */ @@ -556,7 +716,7 @@ void run_mosflm(struct image *image, UnitCell *cell) if ( mosflm->pid == -1 ) { ERROR("Failed to fork for MOSFLM: %s\n", strerror(errno)); free(mosflm); - return; + return 0; } if ( mosflm->pid == 0 ) { @@ -584,7 +744,11 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm->step = 1; /* This starts the "initialisation" procedure */ mosflm->finished_ok = 0; + mosflm->mp = (struct mosflm_private *)ipriv; + mosflm->done = 0; + mosflm->success = 0; + rval = 0; do { fd_set fds; @@ -632,8 +796,50 @@ void run_mosflm(struct image *image, UnitCell *cell) ERROR("MOSFLM doesn't seem to be working properly.\n"); } else { /* Read the mosflm NEWMAT file and get cell if found */ - read_newmat(mosflm->newmatfile, image); + read_newmat(mosflm, mosflm->newmatfile, image); } + rval = mosflm->success; free(mosflm); + return rval; +} + + +IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) +{ + struct mosflm_private *mp; + int need_cell = 0; + + if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1; + if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1; + if ( *indm & INDEXING_USE_LATTICE_TYPE ) need_cell = 1; + + if ( need_cell && (cell == NULL) ) { + ERROR("MOSFLM needs a unit cell for this set of flags.\n"); + return NULL; + } + + /* Flags that MOSFLM knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_CELL_COMBINATIONS + | INDEXING_CHECK_CELL_AXES | INDEXING_CHECK_PEAKS + | INDEXING_USE_LATTICE_TYPE; + + mp = malloc(sizeof(struct mosflm_private)); + if ( mp == NULL ) return NULL; + + mp->ltl = ltl; + mp->template = cell; + mp->indm = *indm; + + return (IndexingPrivate *)mp; +} + + +void mosflm_cleanup(IndexingPrivate *pp) +{ + struct mosflm_private *p; + p = (struct mosflm_private *)pp; + free(p); } diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h index ba1eb73d..a87232b6 100644 --- a/libcrystfel/src/mosflm.h +++ b/libcrystfel/src/mosflm.h @@ -3,13 +3,13 @@ * * Invoke the DPS auto-indexing algorithm through MOSFLM * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: * 2010 Richard Kirian <rkirian@asu.edu> - * 2012 Thomas White <taw@physics.org> + * 2012-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -35,10 +35,16 @@ #include <config.h> #endif -#include "utils.h" +#include "index.h" -extern void run_mosflm(struct image *image, UnitCell *cell); +extern int run_mosflm(struct image *image, IndexingPrivate *ipriv); +extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, + struct detector *det, + struct beam_params *beam, float *ltl); + +extern void mosflm_cleanup(IndexingPrivate *pp); #endif /* MOSFLM_H */ diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index dabc3b1a..6c09053b 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -3,12 +3,12 @@ * * Peak search and other image analysis * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de> * 2011 Andrew Martin <andrew.martin@desy.de> * 2011 Richard Kirian @@ -147,24 +147,15 @@ static int cull_peaks(struct image *image) } -/* cfs, css relative to panel origin */ -static int *make_BgMask(struct image *image, struct panel *p, - double ir_out, double ir_inn) +static void add_crystal_to_mask(struct image *image, struct panel *p, + double ir_inn, int w, int h, + int *mask, Crystal *cr) { Reflection *refl; RefListIterator *iter; - int *mask; - int w, h; - - w = p->max_fs - p->min_fs + 1; - h = p->max_ss - p->min_ss + 1; - mask = calloc(w*h, sizeof(int)); - if ( mask == NULL ) return NULL; - - if ( image->reflections == NULL ) return mask; /* Loop over all reflections */ - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -205,6 +196,27 @@ static int *make_BgMask(struct image *image, struct panel *p, } } +} + + +/* cfs, css relative to panel origin */ +static int *make_BgMask(struct image *image, struct panel *p, double ir_inn) +{ + int *mask; + int w, h; + int i; + + w = p->max_fs - p->min_fs + 1; + h = p->max_ss - p->min_ss + 1; + mask = calloc(w*h, sizeof(int)); + if ( mask == NULL ) return NULL; + + if ( image->crystals == NULL ) return mask; + + for ( i=0; i<image->n_crystals; i++ ) { + add_crystal_to_mask(image, p, ir_inn, + w, h, mask, image->crystals[i]); + } return mask; } @@ -401,9 +413,8 @@ static void search_peaks_in_panel(struct image *image, float threshold, int nrej_dis = 0; int nrej_pro = 0; int nrej_fra = 0; - int nrej_bad = 0; + int nrej_fail = 0; int nrej_snr = 0; - int nrej_sat = 0; int nacc = 0; int ncull; @@ -426,9 +437,6 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Overall threshold */ if ( data[fs+stride*ss] < threshold ) continue; - /* Immediate rejection of pixels above max_adu */ - if ( data[fs+stride*ss] > p->max_adu ) continue; - /* Get gradients */ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss]; dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss]; @@ -494,12 +502,11 @@ static void search_peaks_in_panel(struct image *image, float threshold, /* Centroid peak and get better coordinates. */ r = integrate_peak(image, mask_fs, mask_ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, &saturated); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); - if ( r ) { - /* Bad region - don't detect peak */ - nrej_bad++; - continue; + if ( saturated ) { + image->num_saturated_peaks++; + if ( !use_saturated ) continue; } /* It is possible for the centroid to fall outside the image */ @@ -521,8 +528,9 @@ static void search_peaks_in_panel(struct image *image, float threshold, continue; } - if ( saturated && !use_saturated ) { - nrej_sat++; + if ( r ) { + /* Bad region - don't detect peak */ + nrej_fail++; continue; } @@ -543,10 +551,12 @@ static void search_peaks_in_panel(struct image *image, float threshold, ncull = 0; } + image->num_peaks += nacc; + //STATUS("%i accepted, %i box, %i proximity, %i outside panel, " - // "%i in bad regions, %i with SNR < %g, %i badrow culled, " + // "%i failed integration, %i with SNR < %g, %i badrow culled, " // "%i saturated.\n", - // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad, + // nacc, nrej_dis, nrej_pro, nrej_fra, nrej_fail, // nrej_snr, min_snr, ncull, nrej_sat); if ( ncull != 0 ) { @@ -567,6 +577,8 @@ void search_peaks(struct image *image, float threshold, float min_gradient, image_feature_list_free(image->features); } image->features = image_feature_list_new(); + image->num_peaks = 0; + image->num_saturated_peaks = 0; for ( i=0; i<image->det->n_panels; i++ ) { @@ -581,29 +593,19 @@ void search_peaks(struct image *image, float threshold, float min_gradient, } -double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) +int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst) { - int i; int n_feat = 0; int n_sane = 0; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double min_dist = 0.25; - double stot = 0.0; - - /* Round towards nearest */ - fesetround(1); - - /* Cell basis vectors for this image */ - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + int i; + const double min_dist = 0.25; - /* Loop over peaks, checking proximity to nearest reflection */ for ( i=0; i<image_feature_count(image->features); i++ ) { struct imagefeature *f; struct rvec q; double h,k,l,hd,kd,ld; + int j; /* Assume all image "features" are genuine peaks */ f = image_get_feature(image->features, i); @@ -613,38 +615,42 @@ double peak_lattice_agreement(struct image *image, UnitCell *cell, double *pst) /* Reciprocal space position of found peak */ q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); - /* Decimal and fractional Miller indices of nearest - * reciprocal lattice point */ - hd = q.u * ax + q.v * ay + q.w * az; - kd = q.u * bx + q.v * by + q.w * bz; - ld = q.u * cx + q.v * cy + q.w * cz; - h = lrint(hd); - k = lrint(kd); - l = lrint(ld); - - /* Check distance */ - if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist) - && (fabs(l - ld) < min_dist) ) - { - double sval; - n_sane++; - sval = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0); - stot += 1.0 - sval; - continue; - } + for ( j=0; j<n_cryst; j++ ) { + + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + + cell_get_cartesian(crystal_get_cell(crystals[j]), + &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + + /* Decimal and fractional Miller indices of nearest + * reciprocal lattice point */ + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + /* Check distance */ + if ( (fabs(h - hd) < min_dist) + && (fabs(k - kd) < min_dist) + && (fabs(l - ld) < min_dist) ) + { + n_sane++; + continue; + } - } + } - *pst = stot; - return (double)n_sane / (float)n_feat; -} + } -int peak_sanity_check(struct image *image) -{ - double stot; /* 0 means failed test, 1 means passed test */ - return peak_lattice_agreement(image, image->indexed_cell, &stot) >= 0.5; + return ((double)n_sane / n_feat) >= 0.5; } @@ -705,43 +711,29 @@ static struct integr_ind *sort_reflections(RefList *list, UnitCell *cell, } -/* Integrate the list of predicted reflections in "image" */ -void integrate_reflections(struct image *image, int use_closer, int bgsub, - double min_snr, - double ir_inn, double ir_mid, double ir_out, - int integrate_saturated) +static void integrate_crystal(Crystal *cr, struct image *image, int use_closer, + int bgsub, double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int **bgMasks, + int res_cutoff) { + RefList *reflections; struct integr_ind *il; int n, i; double av = 0.0; int first = 1; - int **bgMasks; + int n_saturated = 0; - if ( num_reflections(image->reflections) == 0 ) return; + reflections = crystal_get_reflections(cr); - il = sort_reflections(image->reflections, image->indexed_cell, &n); + if ( num_reflections(reflections) == 0 ) return; + + il = sort_reflections(reflections, crystal_get_cell(cr), &n); if ( il == NULL ) { ERROR("Couldn't sort reflections\n"); return; } - /* Make background masks for all panels */ - bgMasks = calloc(image->det->n_panels, sizeof(int *)); - if ( bgMasks == NULL ) { - ERROR("Couldn't create list of background masks.\n"); - return; - } - for ( i=0; i<image->det->n_panels; i++ ) { - int *mask; - mask = make_BgMask(image, &image->det->panels[i], - ir_out, ir_inn); - if ( mask == NULL ) { - ERROR("Couldn't create background mask.\n"); - return; - } - bgMasks[i] = mask; - } - for ( i=0; i<n; i++ ) { double fs, ss, intensity; @@ -809,7 +801,7 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, 1, bgMasks[pnum], &saturated); if ( saturated ) { - image->n_saturated++; + n_saturated++; if ( !integrate_saturated ) r = 1; } @@ -840,18 +832,52 @@ void integrate_reflections(struct image *image, int use_closer, int bgsub, } //STATUS("%5.2f A, %5.2f, av %5.2f\n", // 1e10/il[i].res, snr, av); - //if ( av < 1.0 ) break; + if ( res_cutoff && (av < 1.0) ) break; + } + } + + crystal_set_num_saturated_reflections(cr, n_saturated); + crystal_set_resolution_limit(cr, 0.0); + + free(il); +} + + +/* Integrate the list of predicted reflections in "image" */ +void integrate_reflections(struct image *image, int use_closer, int bgsub, + double min_snr, + double ir_inn, double ir_mid, double ir_out, + int integrate_saturated, int res_cutoff) +{ + int i; + int **bgMasks; + + /* Make background masks for all panels */ + bgMasks = calloc(image->det->n_panels, sizeof(int *)); + if ( bgMasks == NULL ) { + ERROR("Couldn't create list of background masks.\n"); + return; + } + for ( i=0; i<image->det->n_panels; i++ ) { + int *mask; + mask = make_BgMask(image, &image->det->panels[i], ir_inn); + if ( mask == NULL ) { + ERROR("Couldn't create background mask.\n"); + return; } + bgMasks[i] = mask; + } + + for ( i=0; i<image->n_crystals; i++ ) { + integrate_crystal(image->crystals[i], image, use_closer, + bgsub, min_snr, ir_inn, ir_mid, ir_out, + integrate_saturated, bgMasks, res_cutoff); } for ( i=0; i<image->det->n_panels; i++ ) { free(bgMasks[i]); } free(bgMasks); - - image->diffracting_resolution = 0.0; - - free(il); } @@ -894,12 +920,17 @@ void validate_peaks(struct image *image, double min_snr, r = integrate_peak(image, f->fs, f->ss, &f_fs, &f_ss, &intensity, &sigma, - ir_inn, ir_mid, ir_out, 0, NULL, &saturated); + ir_inn, ir_mid, ir_out, 1, NULL, &saturated); if ( r ) { n_int++; continue; } + if ( saturated ) { + n_sat++; + if ( !use_saturated ) continue; + } + /* It is possible for the centroid to fall outside the image */ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs) || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) @@ -920,11 +951,6 @@ void validate_peaks(struct image *image, double min_snr, continue; } - if ( saturated && !use_saturated ) { - n_sat++; - continue; - } - /* Add using "better" coordinates */ image_add_feature(flist, f_fs, f_ss, image, intensity, NULL); @@ -936,4 +962,6 @@ void validate_peaks(struct image *image, double min_snr, // n_wtf, n_int, n_dft, n_snr, n_prx, n_sat); image_feature_list_free(image->features); image->features = flist; + image->num_saturated_peaks = n_sat; + image->num_peaks = image_feature_count(flist); } diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h index 3fae13fb..6be728fe 100644 --- a/libcrystfel/src/peaks.h +++ b/libcrystfel/src/peaks.h @@ -3,11 +3,11 @@ * * Peak search and other image analysis * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -45,12 +45,10 @@ extern void search_peaks(struct image *image, float threshold, extern void integrate_reflections(struct image *image, int use_closer, int bgsub, double min_snr, double ir_inn, double ir_mid, double ir_out, - int integrate_saturated); + int integrate_saturated, int res_cutoff); -extern double peak_lattice_agreement(struct image *image, UnitCell *cell, - double *pst); - -extern int peak_sanity_check(struct image *image); +extern int peak_sanity_check(struct image *image, Crystal **crystals, + int n_cryst); extern void estimate_resolution(RefList *list, UnitCell *cell, double *min, double *max); diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c index 7529d12f..3460a4f2 100644 --- a/libcrystfel/src/reax.c +++ b/libcrystfel/src/reax.c @@ -3,11 +3,11 @@ * * A new auto-indexer * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2011-2012 Thomas White <taw@physics.org> + * 2011-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -65,6 +65,9 @@ #define MAX_CANDIDATES (1024) +/* Choose the best solution from this many candidate cells */ +#define MAX_REAX_CELL_CANDIDATES (32) + struct dvec { double x; @@ -103,7 +106,7 @@ struct reax_search struct reax_private { - IndexingPrivate base; + IndexingMethod indm; struct dvec *directions; int n_dir; double angular_inc; @@ -907,7 +910,7 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, } - if ( cl->n_cand >= MAX_CELL_CANDIDATES ) { + if ( cl->n_cand >= MAX_REAX_CELL_CANDIDATES ) { /* "cshift" just fell off the end of the list */ } else { cl->cand[cl->n_cand++] = cshift; @@ -915,15 +918,46 @@ static void add_cell_candidate(struct cell_candidate_list *cl, UnitCell *cnew, } +static int check_cell(struct reax_private *dp, struct image *image, + UnitCell *cell) +{ + UnitCell *out; + Crystal *cr; + + out = cell_new_from_cell(cell); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 0; + } + + crystal_set_cell(cr, out); + + if ( dp->indm & INDEXING_CHECK_PEAKS ) { + if ( !peak_sanity_check(image, &cr, 1) ) { + crystal_free(cr); /* Frees the cell as well */ + return 0; + } + } + + image_add_crystal(image, cr); + + return 1; +} + + static void assemble_cells_from_candidates(struct image *image, struct reax_search *s, - UnitCell *cell) + UnitCell *cell, + struct reax_private *p) { int i, j, k; signed int ti, tj, tk; struct cell_candidate_list cl; - cl.cand = calloc(MAX_CELL_CANDIDATES, sizeof(struct cell_candidate)); + cl.cand = calloc(MAX_REAX_CELL_CANDIDATES, + sizeof(struct cell_candidate)); if ( cl.cand == NULL ) { ERROR("Failed to allocate cell candidate list.\n"); return; @@ -967,7 +1001,9 @@ static void assemble_cells_from_candidates(struct image *image, continue; } - peak_lattice_agreement(image, cnew, &fom); + /* FIXME! */ + //peak_lattice_agreement(image, cnew, &fom); + fom = 1.0; add_cell_candidate(&cl, cnew, fom); } @@ -985,22 +1021,20 @@ static void assemble_cells_from_candidates(struct image *image, cell_get_parameters(cl.cand[i].cell, &a, &b, &c, &al, &be, &ga); cell_get_parameters(cl.cand[i].cell, &aA, &bA, &cA, &alA, &beA, &gaA); - if ( (a - aA) > aA/10.0 ) w = 1; - if ( (b - bA) > bA/10.0 ) w = 1; - if ( (c - cA) > cA/10.0 ) w = 1; - if ( (al - alA) > deg2rad(5.0) ) w = 1; - if ( (be - beA) > deg2rad(5.0) ) w = 1; - if ( (ga - gaA) > deg2rad(5.0) ) w = 1; + if ( fabs(a - aA) > aA/10.0 ) w = 1; + if ( fabs(b - bA) > bA/10.0 ) w = 1; + if ( fabs(c - cA) > cA/10.0 ) w = 1; + if ( fabs(al - alA) > deg2rad(5.0) ) w = 1; + if ( fabs(be - beA) > deg2rad(5.0) ) w = 1; + if ( fabs(ga - gaA) > deg2rad(5.0) ) w = 1; if ( w ) { STATUS("This cell is a long way from that sought:\n"); cell_print(cl.cand[i].cell); } } - image->ncells = cl.n_cand; - assert(image->ncells <= MAX_CELL_CANDIDATES); for ( i=0; i<cl.n_cand; i++ ) { - image->candidate_cells[i] = cl.cand[i].cell; + if ( check_cell(p, image, cl.cand[i].cell) ) break; } free(cl.cand); @@ -1016,7 +1050,6 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) struct reax_search *s; int i; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; fft_in = fftw_malloc(p->nel*sizeof(double)); @@ -1039,7 +1072,7 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) // fft_in, fft_out, p->plan, smin, smax, // image->det, p); - assemble_cells_from_candidates(image, s, cell); + assemble_cells_from_candidates(image, s, cell, p); for ( i=0; i<s->n_search; i++ ) { free(s->search[i].cand); @@ -1051,16 +1084,27 @@ void reax_index(IndexingPrivate *pp, struct image *image, UnitCell *cell) } -IndexingPrivate *reax_prepare() +IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) { struct reax_private *p; int samp; double th; + if ( cell == NULL ) { + ERROR("ReAx needs a unit cell.\n"); + return NULL; + } + p = calloc(1, sizeof(*p)); if ( p == NULL ) return NULL; - p->base.indm = INDEXING_REAX; + /* Flags that ReAx knows about */ + *indm &= INDEXING_METHOD_MASK | INDEXING_CHECK_PEAKS; + + /* Flags that ReAx requires */ + *indm |= INDEXING_USE_LATTICE_TYPE; p->angular_inc = deg2rad(1.0); @@ -1122,6 +1166,8 @@ IndexingPrivate *reax_prepare() p->r_plan = fftw_plan_dft_2d(p->cw, p->ch, p->r_fft_in, p->r_fft_out, 1, FFTW_MEASURE); + p->indm = *indm; + return (IndexingPrivate *)p; } @@ -1130,7 +1176,6 @@ void reax_cleanup(IndexingPrivate *pp) { struct reax_private *p; - assert(pp->indm == INDEXING_REAX); p = (struct reax_private *)pp; free(p->directions); diff --git a/libcrystfel/src/reax.h b/libcrystfel/src/reax.h index 01dfe649..657a2cdf 100644 --- a/libcrystfel/src/reax.h +++ b/libcrystfel/src/reax.h @@ -33,17 +33,26 @@ #include <config.h> #endif +#include "index.h" #include "cell.h" +#include "beam-parameters.h" +#include "detector.h" #ifdef HAVE_FFTW -extern IndexingPrivate *reax_prepare(void); +extern IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl); + extern void reax_cleanup(IndexingPrivate *pp); -extern void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell); + +extern int reax_index(struct image *image, IndexingPrivate *p); #else /* HAVE_FFTW */ -static IndexingPrivate *reax_prepare() +static IndexingPrivate *reax_prepare(IndexingMethod *indm, UnitCell *cell, + const char *filename, struct detector *det, + struct beam_params *beam, float *ltl) { return NULL; } @@ -52,7 +61,7 @@ static void reax_cleanup(IndexingPrivate *pp) { } -static void reax_index(IndexingPrivate *p, struct image *image, UnitCell *cell) +static void reax_index(struct image *image, IndexingPrivate *p) { } diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index b8ff9238..692a860a 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -3,12 +3,12 @@ * * Stream tools * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2011 Andrew Aquila * @@ -37,6 +37,10 @@ #include <stdlib.h> #include <stdio.h> #include <string.h> +#include <assert.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <fcntl.h> #include "cell.h" #include "utils.h" @@ -45,107 +49,17 @@ #include "reflist.h" #include "reflist-utils.h" +#define LATEST_MAJOR_VERSION (2) +#define LATEST_MINOR_VERSION (1) -#define CHUNK_START_MARKER "----- Begin chunk -----" -#define CHUNK_END_MARKER "----- End chunk -----" -#define PEAK_LIST_START_MARKER "Peaks from peak search" -#define PEAK_LIST_END_MARKER "End of peak list" -#define REFLECTION_START_MARKER "Reflections measured after indexing" -/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also - * used to terminate a standalone list of reflections */ -static void exclusive(const char *a, const char *b) +struct _stream { - ERROR("The stream options '%s' and '%s' are mutually exclusive.\n", - a, b); -} - - -int parse_stream_flags(const char *a) -{ - int n, i; - int ret = STREAM_NONE; - char **flags; - - n = assplode(a, ",", &flags, ASSPLODE_NONE); - - for ( i=0; i<n; i++ ) { - - if ( strcmp(flags[i], "integrated") == 0) { - - ret |= STREAM_INTEGRATED; - - } else if ( strcmp(flags[i], "peaks") == 0) { - if ( ret & STREAM_PEAKS_IF_INDEXED ) { - exclusive("peaks", "peaksifindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) { - exclusive("peaks", "peaksifnotindexed"); - return -1; - } - ret |= STREAM_PEAKS; - - } else if ( strcmp(flags[i], "peaksifindexed") == 0) { - if ( ret & STREAM_PEAKS ) { - exclusive("peaks", "peaksifindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) { - exclusive("peaksifnotindexed", - "peaksifindexed"); - return -1; - } - ret |= STREAM_PEAKS_IF_INDEXED; - - } else if ( strcmp(flags[i], "peaksifnotindexed") == 0) { - if ( ret & STREAM_PEAKS ) { - exclusive("peaks", "peaksifnotindexed"); - return -1; - } - if ( ret & STREAM_PEAKS_IF_INDEXED ) { - exclusive("peaksifnotindexed", - "peaksifindexed"); - return -1; - } - ret |= STREAM_PEAKS_IF_NOT_INDEXED; - - } else { - ERROR("Unrecognised stream flag '%s'\n", flags[i]); - return -1; - } - - free(flags[i]); - - } - free(flags); - - return ret; -} - - -int count_patterns(FILE *fh) -{ - char *rval; - - int n_total_patterns = 0; - do { - char line[1024]; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++; - - } while ( rval != NULL ); - - if ( ferror(fh) ) { - ERROR("Read error while counting patterns.\n"); - return 0; - } + FILE *fh; - return n_total_patterns; -} + int major_version; + int minor_version; +}; static int read_peaks(FILE *fh, struct image *image) @@ -215,96 +129,107 @@ static void write_peaks(struct image *image, FILE *ofh) } -void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f) +static void write_crystal(Stream *st, Crystal *cr, int include_reflections) { + UnitCell *cell; + RefList *reflist; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double a, b, c, al, be, ga; - fprintf(ofh, CHUNK_START_MARKER"\n"); + fprintf(st->fh, CRYSTAL_START_MARKER"\n"); - fprintf(ofh, "Image filename: %s\n", i->filename); + cell = crystal_get_cell(cr); + assert(cell != NULL); - if ( i->indexed_cell != NULL ) { + cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga); + fprintf(st->fh, "Cell parameters %7.5f %7.5f %7.5f nm," + " %7.5f %7.5f %7.5f deg\n", + a*1.0e9, b*1.0e9, c*1.0e9, + rad2deg(al), rad2deg(be), rad2deg(ga)); - cell_get_parameters(i->indexed_cell, &a, &b, &c, - &al, &be, &ga); - fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm," - " %7.5f %7.5f %7.5f deg\n", - a*1.0e9, b*1.0e9, c*1.0e9, - rad2deg(al), rad2deg(be), rad2deg(ga)); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + fprintf(st->fh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", + asx/1e9, asy/1e9, asz/1e9); + fprintf(st->fh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + bsx/1e9, bsy/1e9, bsz/1e9); + fprintf(st->fh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", + csx/1e9, csy/1e9, csz/1e9); - cell_get_reciprocal(i->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n", - asx/1e9, asy/1e9, asz/1e9); - fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - bsx/1e9, bsy/1e9, bsz/1e9); - fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n", - csx/1e9, csy/1e9, csz/1e9); + reflist = crystal_get_reflections(cr); + if ( reflist != NULL ) { - } else { + fprintf(st->fh, "diffraction_resolution_limit" + " = %.2f nm^-1 or %.2f A\n", + crystal_get_resolution_limit(cr)/1e9, + 1e9 / crystal_get_resolution_limit(cr)); - fprintf(ofh, "No unit cell from indexing.\n"); + fprintf(st->fh, "num_saturated_reflections = %lli\n", + crystal_get_num_saturated_reflections(cr)); } - fprintf(ofh, "photon_energy_eV = %f\n", - J_to_eV(ph_lambda_to_en(i->lambda))); + if ( include_reflections ) { + + if ( reflist != NULL ) { - if ( i->reflections != NULL ) { + fprintf(st->fh, REFLECTION_START_MARKER"\n"); + write_reflections_to_file(st->fh, reflist); + fprintf(st->fh, REFLECTION_END_MARKER"\n"); - fprintf(ofh, "diffraction_resolution_limit" - " = %.2f nm^-1 or %.2f A\n", - i->diffracting_resolution/1e9, - 1e9 / i->diffracting_resolution); + } else { - fprintf(ofh, "num_saturated_reflections" - " = %lli\n", i->n_saturated); + fprintf(st->fh, "No integrated reflections.\n"); + } } + fprintf(st->fh, CRYSTAL_END_MARKER"\n"); +} + + +void write_chunk(Stream *st, struct image *i, struct hdfile *hdfile, + int include_peaks, int include_reflections) +{ + int j; + + fprintf(st->fh, CHUNK_START_MARKER"\n"); + + fprintf(st->fh, "Image filename: %s\n", i->filename); + fprintf(st->fh, "indexed_by = %s\n", indexer_str(i->indexed_by)); + if ( i->det != NULL ) { int j; for ( j=0; j<i->det->n_panels; j++ ) { - fprintf(ofh, "camera_length_%s = %f\n", + fprintf(st->fh, "camera_length_%s = %f\n", i->det->panels[j].name, i->det->panels[j].clen); } } - copy_hdf5_fields(hdfile, i->copyme, ofh); + copy_hdf5_fields(hdfile, i->copyme, st->fh); - if ( (f & STREAM_PEAKS) - || ((f & STREAM_PEAKS_IF_INDEXED) && (i->indexed_cell != NULL)) - || ((f & STREAM_PEAKS_IF_NOT_INDEXED) && (i->indexed_cell == NULL)) ) - { - fprintf(ofh, "\n"); - write_peaks(i, ofh); + fprintf(st->fh, "num_peaks = %lli\n", i->num_peaks); + fprintf(st->fh, "num_saturated_peaks = %lli\n", i->num_saturated_peaks); + if ( include_peaks ) { + write_peaks(i, st->fh); } - if ( f & STREAM_INTEGRATED ) { - - fprintf(ofh, "\n"); - - if ( i->reflections != NULL ) { - - fprintf(ofh, REFLECTION_START_MARKER"\n"); - write_reflections_to_file(ofh, i->reflections); - fprintf(ofh, REFLECTION_END_MARKER"\n"); - - } else { - - fprintf(ofh, "No integrated reflections.\n"); + fprintf(st->fh, "photon_energy_eV = %f\n", + J_to_eV(ph_lambda_to_en(i->lambda))); - } + for ( j=0; j<i->n_crystals; j++ ) { + write_crystal(st, i->crystals[j], include_reflections); } - fprintf(ofh, CHUNK_END_MARKER"\n\n"); + fprintf(st->fh, CHUNK_END_MARKER"\n"); + + fflush(st->fh); } @@ -328,8 +253,7 @@ static int find_start_of_chunk(FILE *fh) } -/* Read the next chunk from a stream and fill in 'image' */ -int read_chunk(FILE *fh, struct image *image) +void read_crystal(Stream *st, struct image *image) { char line[1024]; char *rval = NULL; @@ -337,22 +261,116 @@ int read_chunk(FILE *fh, struct image *image) int have_as = 0; int have_bs = 0; int have_cs = 0; + Crystal *cr; + int n; + Crystal **crystals_new; + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal!\n"); + return; + } + + do { + + float u, v, w; + + rval = fgets(line, 1023, st->fh); + + /* Trouble? */ + if ( rval == NULL ) break; + + chomp(line); + if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) { + as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; + have_as = 1; + } + + if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) { + bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; + have_bs = 1; + } + + if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) { + cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; + have_cs = 1; + } + + if ( have_as && have_bs && have_cs ) { + + UnitCell *cell; + + cell = crystal_get_cell(cr); + + if ( cell != NULL ) { + ERROR("Duplicate cell found in stream!\n"); + ERROR("I'll use the most recent one.\n"); + cell_free(cell); + } + + cell = cell_new_from_reciprocal_axes(as, bs, cs); + crystal_set_cell(cr, cell); + + have_as = 0; have_bs = 0; have_cs = 0; + + } + + if ( strncmp(line, "num_saturated_reflections = ", 28) == 0 ) { + int n = atoi(line+28); + crystal_set_num_saturated_reflections(cr, n); + } + + if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { + + RefList *reflections; + + reflections = read_reflections_from_file(st->fh); + + if ( reflections == NULL ) { + ERROR("Failed while reading reflections\n"); + break; + } + + crystal_set_reflections(cr, reflections); + + } + + if ( strcmp(line, CRYSTAL_END_MARKER) == 0 ) break; + + } while ( 1 ); + + /* Add crystal to the list for this image */ + n = image->n_crystals+1; + crystals_new = realloc(image->crystals, n*sizeof(Crystal *)); + + if ( crystals_new == NULL ) { + ERROR("Failed to expand crystal list!\n"); + } else { + image->crystals = crystals_new; + image->crystals[image->n_crystals++] = cr; + } + +} + + +/* Read the next chunk from a stream and fill in 'image' */ +int read_chunk(Stream *st, struct image *image) +{ + char line[1024]; + char *rval = NULL; int have_filename = 0; int have_ev = 0; - if ( find_start_of_chunk(fh) ) return 1; + if ( find_start_of_chunk(st->fh) ) return 1; image->lambda = -1.0; image->features = NULL; - image->reflections = NULL; - image->indexed_cell = NULL; - image->n_saturated = 0; + image->crystals = NULL; + image->n_crystals = 0; do { - float u, v, w; - - rval = fgets(line, 1023, fh); + rval = fgets(line, 1023, st->fh); /* Trouble? */ if ( rval == NULL ) break; @@ -364,6 +382,11 @@ int read_chunk(FILE *fh, struct image *image) have_filename = 1; } + if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) { + image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19))); + have_ev = 1; + } + if ( strncmp(line, "camera_length_", 14) == 0 ) { if ( image->det != NULL ) { @@ -390,56 +413,19 @@ int read_chunk(FILE *fh, struct image *image) } } - if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) { - as.u = u*1e9; as.v = v*1e9; as.w = w*1e9; - have_as = 1; - } - - if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) { - bs.u = u*1e9; bs.v = v*1e9; bs.w = w*1e9; - have_bs = 1; - } - - if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) { - cs.u = u*1e9; cs.v = v*1e9; cs.w = w*1e9; - have_cs = 1; - } - - if ( have_as && have_bs && have_cs ) { - if ( image->indexed_cell != NULL ) { - ERROR("Duplicate cell found in stream!\n"); - cell_free(image->indexed_cell); - } - image->indexed_cell = cell_new_from_reciprocal_axes(as, - bs, - cs); - have_as = 0; have_bs = 0; have_cs = 0; - } - - if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) { - image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19))); - have_ev = 1; - } - - if ( strncmp(line, "num_saturated_reflections = ", 28) == 0 ) { - image->n_saturated = atoi(line+28); - } - if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) { - if ( read_peaks(fh, image) ) { + if ( read_peaks(st->fh, image) ) { ERROR("Failed while reading peaks\n"); return 1; } } - if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { - image->reflections = read_reflections_from_file(fh); - if ( image->reflections == NULL ) { - ERROR("Failed while reading reflections\n"); - return 1; - } + if ( strcmp(line, CRYSTAL_START_MARKER) == 0 ) { + read_crystal(st, image); } + /* A chunk must have at least a filename and a wavelength, + * otherwise it's incomplete */ if ( strcmp(line, CHUNK_END_MARKER) == 0 ) { if ( have_filename && have_ev ) return 0; ERROR("Incomplete chunk found in input file.\n"); @@ -448,6 +434,10 @@ int read_chunk(FILE *fh, struct image *image) } while ( 1 ); + if ( !feof(st->fh) ) { + ERROR("Error reading stream.\n"); + } + return 1; /* Either error or EOF, don't care because we will complain * on the terminal if it was an error. */ } @@ -457,7 +447,6 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) { int i; - fprintf(ofh, "CrystFEL stream format 2.0\n"); fprintf(ofh, "Command line:"); for ( i=0; i<argc; i++ ) { fprintf(ofh, " %s", argv[i]); @@ -467,25 +456,87 @@ void write_stream_header(FILE *ofh, int argc, char *argv[]) } -int skip_some_files(FILE *fh, int n) +Stream *open_stream_for_read(const char *filename) { - char *rval = NULL; - int n_patterns = 0; + Stream *st; - do { + st = malloc(sizeof(struct _stream)); + if ( st == NULL ) return NULL; - char line[1024]; + st->fh = fopen(filename, "r"); + if ( st->fh == NULL ) { + free(st); + return NULL; + } - if ( n_patterns == n ) return 0; + char line[1024]; + char *rval; - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++; + rval = fgets(line, 1023, st->fh); + if ( rval == NULL ) { + ERROR("Failed to read stream version.\n"); + close_stream(st); + return NULL; + } - } while ( rval != NULL ); + if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) { + st->major_version = 2; + st->minor_version = 0; + } else if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) { + st->major_version = 2; + st->minor_version = 1; + } else { + ERROR("Invalid stream, or stream format is too new.\n"); + close_stream(st); + return NULL; + } - return 1; + return st; +} + + +Stream *open_stream_fd_for_write(int fd) +{ + Stream *st; + + st = malloc(sizeof(struct _stream)); + if ( st == NULL ) return NULL; + + st->fh = fdopen(fd, "w"); + if ( st->fh == NULL ) { + free(st); + return NULL; + } + + st->major_version = LATEST_MAJOR_VERSION; + st->minor_version = LATEST_MINOR_VERSION; + + fprintf(st->fh, "CrystFEL stream format %i.%i\n", + st->major_version, st->minor_version); + + return st; +} + + + +Stream *open_stream_for_write(const char *filename) +{ + int fd; + + fd = open(filename, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR); + if ( fd == -1 ) { + ERROR("Failed to open stream.\n"); + return NULL; + } + + return open_stream_fd_for_write(fd); +} + + +void close_stream(Stream *st) +{ + fclose(st->fh); + free(st); } @@ -493,20 +544,53 @@ int is_stream(const char *filename) { FILE *fh; char line[1024]; - char *rval = NULL; + char *rval; + fh = fopen(filename, "r"); - if ( fh == NULL ) { - return -1; - } + if ( fh == NULL ) return 0; + rval = fgets(line, 1023, fh); fclose(fh); - if ( rval == NULL ) { - return -1; - } - if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) { - return 1; - } - else { - return 0; + if ( rval == NULL ) return 0; + + if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) return 1; + if ( strncmp(line, "CrystFEL stream format 2.1", 26) == 0 ) return 1; + + return 0; +} + + +void write_line(Stream *st, const char *line) +{ + fprintf(st->fh, "%s\n", line); +} + + +void write_command(Stream *st, int argc, char *argv[]) +{ + int i; + + for ( i=0; i<argc; i++ ) { + if ( i > 0 ) fprintf(st->fh, " "); + fprintf(st->fh, "%s", argv[i]); } + fprintf(st->fh, "\n"); +} + + +/** + * rewind_stream: + * @st: A %Stream + * + * Attempts to set the file pointer for @st to the start of the stream, so that + * later calls to read_chunk() will repeat the sequence of chunks from the + * start. + * + * Programs must not assume that this operation always succeeds! + * + * Returns: non-zero if the stream could not be rewound. + */ +int rewind_stream(Stream *st) +{ + return fseek(st->fh, 0, SEEK_SET); } diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index a544798f..08e680bc 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -3,11 +3,11 @@ * * Stream tools * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Andrew Aquila * * This file is part of CrystFEL. @@ -38,30 +38,31 @@ struct image; struct hdfile; -/* Possible options dictating what goes into the output stream */ -enum -{ - STREAM_NONE = 0, - STREAM_INTEGRATED = 1<<0, - STREAM_PEAKS = 1<<2, - STREAM_PEAKS_IF_INDEXED = 1<<3, - STREAM_PEAKS_IF_NOT_INDEXED = 1<<4, -}; +#define CHUNK_START_MARKER "----- Begin chunk -----" +#define CHUNK_END_MARKER "----- End chunk -----" +#define PEAK_LIST_START_MARKER "Peaks from peak search" +#define PEAK_LIST_END_MARKER "End of peak list" +#define CRYSTAL_START_MARKER "--- Begin crystal" +#define CRYSTAL_END_MARKER "--- End crystal" +#define REFLECTION_START_MARKER "Reflections measured after indexing" +/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also + * used to terminate a standalone list of reflections */ +typedef struct _stream Stream; -extern int count_patterns(FILE *fh); +extern Stream *open_stream_for_read(const char *filename); +extern Stream *open_stream_for_write(const char *filename); +extern Stream *open_stream_fd_for_write(int fd); +extern void close_stream(Stream *st); -extern void write_stream_header(FILE *ofh, int argc, char *argv[]); +extern int read_chunk(Stream *st, struct image *image); +extern void write_chunk(Stream *st, struct image *image, struct hdfile *hdfile, + int include_peaks, int include_reflections); -extern void write_chunk(FILE *ofh, struct image *image, struct hdfile *hdfile, - int flags); - -extern int parse_stream_flags(const char *a); - -extern int read_chunk(FILE *fh, struct image *image); - -extern int skip_some_files(FILE *fh, int n); +extern void write_line(Stream *st, const char *line); +extern void write_command(Stream *st, int argc, char *argv[]); +extern int rewind_stream(Stream *st); extern int is_stream(const char *filename); #endif /* STREAM_H */ diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index a652e413..3c2011d1 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -1069,11 +1069,12 @@ static void do_op(const IntegerMatrix *op, * This function applies the @idx-th symmetry operation from @ops to the * reflection @h, @k, @l, and stores the result at @he, @ke and @le. * - * If you don't mind that the same equivalent might appear twice, simply call - * this function the number of times returned by num_ops(), using the actual - * point group. If repeating the same equivalent twice (for example, if the - * given reflection is a special high-symmetry one), call special_position() - * first to get a "specialised" SymOpList and use that instead. + * Call this function multiple times with idx=0 .. num_equivs(ops, m) to get all + * of the equivalent reflections in turn. + * + * If you don't mind that the same equivalent might appear twice, simply let + * @m = NULL. Otherwise, call new_symopmask() and then special_position() to + * set up a %SymOpMask appropriately. **/ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx, signed int h, signed int k, signed int l, @@ -1133,8 +1134,8 @@ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx, * @k: index of a reflection * @l: index of a reflection * - * This function determines which operations in @ops map the reflection @h, @k, - * @l onto itself, and uses @m to flag the operations in @ops which cause this. + * This function sets up @m to contain information about which operations in + * @ops map the reflection @h, @k, @l onto itself. * **/ void special_position(const SymOpList *ops, SymOpMask *m, @@ -1467,7 +1468,7 @@ static SymOpList *flack_reorder(const SymOpList *source) * operators, which are the symmetry operations which can be added to @target * to generate @source. Only rotations are allowable - no mirrors nor * inversions. - * To count the number of possibilities, use num_ops() on the result. + * To count the number of possibilities, use num_equivs() on the result. * * The algorithm used is "Algorithm A" from Flack (1987), Acta Cryst A43 p564. * diff --git a/scripts/alternate-stream b/scripts/alternate-stream index 1ac2fd3a..044edb24 100755 --- a/scripts/alternate-stream +++ b/scripts/alternate-stream @@ -9,6 +9,10 @@ open(OFH_TWO, "> ".$ARGV[2]); my $line; my $alt = 0; +my $header = <FH>; +printf(OFH_ONE $header); +printf(OFH_TWO $header); + while ( $line = <FH> ) { if ( $line =~ /^-----\ Begin chunk\ -----$/ ) { diff --git a/scripts/create-xscale b/scripts/create-xscale index 17dc6280..5d810f07 100755 --- a/scripts/create-xscale +++ b/scripts/create-xscale @@ -21,7 +21,7 @@ while ( $line = <FH> ) { chomp($line); - if ( $line =~ /^\s+([0-9]+)\s+([0-9]+)\s+([0-9]+)\s+([0-9\.\-]+)/ ) { + if ( $line =~ /^\s+([0-9]+\-)\s+([0-9]+\-)\s+([0-9]+\-)\s+([0-9\.\-]+)/ ) { my $h = $1; my $k = $2; diff --git a/scripts/peak-intensity b/scripts/peak-intensity index 46f5be9e..13bed699 100755 --- a/scripts/peak-intensity +++ b/scripts/peak-intensity @@ -8,6 +8,9 @@ my $line; my $total_i; my $n = 0; my $n_patt = 0; +my $num_peaks_tot = 0; +my $num_sat_peaks_tot = 0; +my $num_pats_with_sat = 0; while ( $line = <FH> ) { @@ -27,9 +30,21 @@ while ( $line = <FH> ) { } + if ( $line =~ /^num_saturated_peaks\s=\s(\d+)$/ ) { + $num_sat_peaks_tot += $1; + if ( $1 > 0 ) { + $num_pats_with_sat++; + } + } + } printf("%i patterns, %i peaks, %.2f total ADU\n", $n_patt, $n, $total_i); printf("Mean %i peaks per hit\n", $n / $n_patt); printf("Mean %.2f ADU per peak\n", $total_i / $n); printf("Mean %.2f ADU per hit\n", $total_i / $n_patt); +printf("%i out of %i patterns contained any saturation\n", $num_pats_with_sat, $n_patt); +if ( $num_pats_with_sat > 0 ) { + printf(" of those, there was an average of %.2f saturated peaks per pattern\n", + $num_sat_peaks_tot/$num_pats_with_sat); +} diff --git a/src/check_hkl.c b/src/check_hkl.c index 8e171904..f82855ab 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -380,6 +380,9 @@ int main(int argc, char *argv[]) } break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; diff --git a/src/compare_hkl.c b/src/compare_hkl.c index f7776598..c3e7a547 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -95,6 +95,7 @@ static void show_help(const char *s) " --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n" " --rmin=<res> Set a lower resolution limit (m^-1).\n" " --rmax=<res> Set an upper resolution limit (m^-1).\n" +" --intensity-shells Use shells of intensity instead of resolution.\n" "\n" " -h, --help Display this help message.\n" ); @@ -770,6 +771,9 @@ int main(int argc, char *argv[]) shell_file = strdup(optarg); break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c index 19aa6f35..ee3b55f0 100644 --- a/src/dw-hdfsee.c +++ b/src/dw-hdfsee.c @@ -1819,7 +1819,6 @@ DisplayWindow *displaywindow_open(const char *filename, const char *peaks, dw->hdfile = hdfile_open(filename); if ( dw->hdfile == NULL ) { - ERROR("Couldn't open file '%s'\n", filename); free(dw); return NULL; } else { diff --git a/src/get_hkl.c b/src/get_hkl.c index 7b66e427..c879a51a 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -425,6 +425,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; diff --git a/src/hdfsee.c b/src/hdfsee.c index e5d90456..894b8af9 100644 --- a/src/hdfsee.c +++ b/src/hdfsee.c @@ -223,6 +223,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; @@ -263,11 +266,7 @@ int main(int argc, char *argv[]) ring_radii, n_rings, ring_size); - if ( main_window_list[i] == NULL ) { - ERROR("Couldn't open display window\n"); - } else { - main_n_windows++; - } + if ( main_window_list[i] != NULL ) main_n_windows++; } if ( main_n_windows == 0 ) return 0; diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index f86ae0bb..0498a532 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -57,7 +57,7 @@ struct scale_queue_args { RefList *reference; - struct image *images; + Crystal **crystals; int n_started; double max_shift; }; @@ -65,7 +65,7 @@ struct scale_queue_args struct scale_worker_args { - struct image *image; + Crystal *crystal; double shift; RefList *reference; }; @@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs) wargs = malloc(sizeof(struct scale_worker_args)); wargs->reference = qargs->reference; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs) static void run_scale_job(void *vwargs, int cookie) { struct scale_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *reference = wargs->reference; Reflection *refl; RefListIterator *iter; double num = 0.0; double den = 0.0; double corr; + const double osf = crystal_get_osf(cr); - if ( image->pr_dud ) { + if ( crystal_get_user_flag(cr) ) { wargs->shift = 0.0; return; } - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie) Ihl = get_intensity(refl) / get_partiality(refl); esd = get_esd_intensity(refl) / get_partiality(refl); - num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0); - den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0); + num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); + den += pow(Ih, 2.0)/pow(esd/osf, 2.0); } @@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie) corr = num / den; if ( !isnan(corr) && !isinf(corr) ) { - image->osf *= corr; + crystal_set_osf(cr, osf*corr); } wargs->shift = fabs(corr-1.0); @@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs) } -static double iterate_scale(struct image *images, int n, RefList *reference, +static double iterate_scale(Crystal **crystals, int n, RefList *reference, int n_threads) { struct scale_queue_args qargs; @@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference, qargs.reference = reference; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; qargs.max_shift = 0.0; run_threads(n_threads, run_scale_job, create_scale_job, @@ -178,14 +179,14 @@ struct merge_queue_args { RefList *full; pthread_mutex_t full_lock; - struct image *images; + Crystal **crystals; int n_started; }; struct merge_worker_args { - struct image *image; + Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; }; @@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs) wargs->full = qargs->full; wargs->full_lock = &qargs->full_lock; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs) static void run_merge_job(void *vwargs, int cookie) { struct merge_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr)) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(struct image *images, int n, int n_threads) +static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) { RefList *full; struct merge_queue_args qargs; @@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; pthread_mutex_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) struct esd_queue_args { RefList *full; - struct image *images; + Crystal **crystals; int n_started; }; struct esd_worker_args { - struct image *image; + Crystal *crystal; RefList *full; }; @@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs) wargs = malloc(sizeof(struct esd_worker_args)); wargs->full = qargs->full; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs) static void run_esd_job(void *vwargs, int cookie) { struct esd_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr) ) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs) } -static void calculate_esds(struct image *images, int n, RefList *full, +static void calculate_esds(Crystal **crystals, int n, RefList *full, int n_threads, int min_red) { struct esd_queue_args qargs; @@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; for ( refl = first_refl(full, &iter); refl != NULL; @@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, /* Scale the stack of images */ -RefList *scale_intensities(struct image *images, int n, RefList *gref, +RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, int n_threads, int noscale) { int i; @@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, RefList *full = NULL; const int min_redundancy = 3; - for ( i=0; i<n; i++ ) images[i].osf = 1.0; + for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0); if ( noscale ) { - full = lsq_intensities(images, n, n_threads); - calculate_esds(images, n, full, n_threads, min_redundancy); + full = lsq_intensities(crystals, n, n_threads); + calculate_esds(crystals, n, full, n_threads, min_redundancy); return full; } /* No reference -> create an initial list to refine against */ if ( gref == NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* Iterate */ @@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, reference = full; } - max_corr = iterate_scale(images, n, reference, n_threads); + max_corr = iterate_scale(crystals, n, reference, n_threads); //STATUS("Scaling iteration %2i: max correction = %5.2f\n", // i+1, max_corr); /* No reference -> generate list for next iteration */ if ( gref == NULL ) { reflist_free(full); - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } //show_scale_factors(images, n); @@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); if ( gref != NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* else we already did it */ - calculate_esds(images, n, full, n_threads, min_redundancy); + calculate_esds(crystals, n, full, n_threads, min_redundancy); return full; } diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h index 5425b1ff..80940347 100644 --- a/src/hrs-scaling.h +++ b/src/hrs-scaling.h @@ -35,9 +35,10 @@ #endif -#include "image.h" +#include "crystal.h" +#include "reflist.h" -extern RefList *scale_intensities(struct image *images, int n, +extern RefList *scale_intensities(Crystal **crystals, int n, RefList *reference, int n_threads, int noscale); diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 32bcba15..8da42d32 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -73,14 +73,27 @@ #define STATS_EVERY_N_SECONDS (5) +/* Information about the indexing process for one pattern */ +struct pattern_args +{ + /* "Input" */ + char *filename; + + /* "Output" */ + int n_crystals; +}; + + struct sandbox { pthread_mutex_t lock; - int n_indexable; int n_processed; - int n_indexable_last_stats; + int n_hadcrystals; + int n_crystals; int n_processed_last_stats; + int n_hadcrystals_last_stats; + int n_crystals_last_stats; int t_last_stats; struct index_args *iargs; @@ -91,6 +104,7 @@ struct sandbox FILE **fhs; int *running; + int *waiting; FILE **result_fhs; int *filename_pipes; int *stream_pipe_read; @@ -164,31 +178,26 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, static void process_image(const struct index_args *iargs, - struct pattern_args *pargs, FILE *ofh, + struct pattern_args *pargs, Stream *st, int cookie) { float *data_for_measurement; size_t data_size; - UnitCell *cell = iargs->cell; - int config_cmfilter = iargs->config_cmfilter; - int config_noisefilter = iargs->config_noisefilter; - int config_verbose = iargs->config_verbose; - IndexingMethod *indm = iargs->indm; int check; struct hdfile *hdfile; struct image image; + int i; image.features = NULL; image.data = NULL; image.flags = NULL; - image.indexed_cell = NULL; image.copyme = iargs->copyme; - image.reflections = NULL; - image.n_saturated = 0; image.id = cookie; image.filename = pargs->filename; image.beam = iargs->beam; image.det = iargs->det; + image.crystals = NULL; + image.n_crystals = 0; hdfile = hdfile_open(image.filename); if ( hdfile == NULL ) return; @@ -247,16 +256,14 @@ static void process_image(const struct index_args *iargs, return; } - if ( config_cmfilter ) { - filter_cm(&image); - } + if ( iargs->cmfilter ) filter_cm(&image); /* Take snapshot of image after CM subtraction but before * the aggressive noise filter. */ data_size = image.width * image.height * sizeof(float); data_for_measurement = malloc(data_size); - if ( config_noisefilter ) { + if ( iargs->noisefilter ) { filter_noise(&image, data_for_measurement); } else { memcpy(data_for_measurement, image.data, data_size); @@ -291,41 +298,54 @@ static void process_image(const struct index_args *iargs, free(image.data); image.data = data_for_measurement; - /* Calculate orientation matrix (by magic) */ + /* Index the pattern */ + index_pattern(&image, iargs->indm, iargs->ipriv); + + pargs->n_crystals = image.n_crystals; + + /* Default beam parameters */ image.div = image.beam->divergence; image.bw = image.beam->bandwidth; - image.profile_radius = image.beam->profile_radius; - - index_pattern(&image, cell, indm, iargs->cellr, - config_verbose, iargs->ipriv, - iargs->config_insane, iargs->tols); - - if ( image.indexed_cell != NULL ) { - pargs->indexable = 1; - image.reflections = find_intersections(&image, - image.indexed_cell); - if (image.reflections != NULL) { - integrate_reflections(&image, - iargs->config_closer, - iargs->config_bgsub, - iargs->min_int_snr, - iargs->ir_inn, - iargs->ir_mid, - iargs->ir_out, - iargs->integrate_saturated); + + /* Integrate each crystal's diffraction spots */ + for ( i=0; i<image.n_crystals; i++ ) { + + RefList *reflections; + + /* Set default crystal parameter(s) */ + crystal_set_profile_radius(image.crystals[i], + image.beam->profile_radius); + + if ( iargs->integrate_found ) { + reflections = select_intersections(&image, + image.crystals[i]); + } else { + reflections = find_intersections(&image, + image.crystals[i]); } - } else { - image.reflections = NULL; - } - write_chunk(ofh, &image, hdfile, iargs->stream_flags); - fprintf(ofh, "END\n"); - fflush(ofh); + crystal_set_reflections(image.crystals[i], reflections); + + } - /* Only free cell if found */ - cell_free(image.indexed_cell); + /* Integrate all the crystals at once - need all the crystals so that + * overlaps can be detected. */ + integrate_reflections(&image, iargs->closer, + iargs->bgsub, + iargs->min_int_snr, + iargs->ir_inn, + iargs->ir_mid, + iargs->ir_out, + iargs->integrate_saturated, + iargs->res_cutoff); + + write_chunk(st, &image, hdfile, + iargs->stream_peaks, iargs->stream_refls); + + for ( i=0; i<image.n_crystals; i++ ) { + crystal_free(image.crystals[i]); + } - reflist_free(image.reflections); free(image.data); if ( image.flags != NULL ) free(image.flags); image_feature_list_free(image.features); @@ -334,7 +354,7 @@ static void process_image(const struct index_args *iargs, static void run_work(const struct index_args *iargs, - int filename_pipe, int results_pipe, FILE *ofh, + int filename_pipe, int results_pipe, Stream *st, int cookie) { int allDone = 0; @@ -380,12 +400,12 @@ static void run_work(const struct index_args *iargs, } else { pargs.filename = line; - pargs.indexable = 0; + pargs.n_crystals = 0; - process_image(iargs, &pargs, ofh, cookie); + process_image(iargs, &pargs, st, cookie); /* Request another image */ - c = sprintf(buf, "%i\n", pargs.indexable); + c = sprintf(buf, "%i\n", pargs.n_crystals); w = write(results_pipe, buf, c); if ( w < 0 ) { ERROR("write P0\n"); @@ -397,7 +417,9 @@ static void run_work(const struct index_args *iargs, } - cleanup_indexing(iargs->ipriv); + write_line(st, "DONE"); + + cleanup_indexing(iargs->indm, iargs->ipriv); free(iargs->indm); free(iargs->ipriv); free_detector_geometry(iargs->det); @@ -462,15 +484,26 @@ static int pump_chunk(FILE *fh, FILE *ofh) } - if ( strcmp(line, "END\n") == 0 ) { + if ( strcmp(line, "FLUSH\n") == 0 ) { chunk_finished = 1; - } else { + continue; + } + + if ( strcmp(line, "DONE\n") == 0 ) { + return 1; + } + + fprintf(ofh, "%s", line); + + if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) { + chunk_finished = 1; + } + if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) { chunk_started = 1; - fprintf(ofh, "%s", line); } - } while ( !chunk_finished ); + } while ( !chunk_finished ); return 0; } @@ -524,10 +557,13 @@ static void *run_reader(void *sbv) if ( !sb->running[i] ) continue; - if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) continue; + if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) { + continue; + } if ( pump_chunk(sb->fhs[i], sb->ofh) ) { sb->running[i] = 0; + sb->waiting[i] = 0; } } @@ -544,7 +580,8 @@ static void *run_reader(void *sbv) } -static void start_worker_process(struct sandbox *sb, int slot) +static void start_worker_process(struct sandbox *sb, int slot, + int argc, char *argv[]) { pid_t p; int filename_pipe[2]; @@ -568,7 +605,7 @@ static void start_worker_process(struct sandbox *sb, int slot) if ( p == 0 ) { - FILE *sfh; + Stream *st; int j; struct sigaction sa; int r; @@ -605,10 +642,12 @@ static void start_worker_process(struct sandbox *sb, int slot) close(filename_pipe[1]); close(result_pipe[0]); - sfh = fdopen(sb->stream_pipe_write[slot], "w"); + st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); + write_command(st, argc, argv); + write_line(st, "FLUSH"); run_work(sb->iargs, filename_pipe[0], result_pipe[1], - sfh, slot); - fclose(sfh); + st, slot); + close_stream(st); //close(filename_pipe[0]); close(result_pipe[1]); @@ -654,6 +693,7 @@ static void handle_zombie(struct sandbox *sb) int status, p; if ( !sb->running[i] ) continue; + if ( sb->waiting[i] ) continue; p = waitpid(sb->pids[i], &status, WNOHANG); @@ -665,6 +705,7 @@ static void handle_zombie(struct sandbox *sb) if ( p == sb->pids[i] ) { sb->running[i] = 0; + sb->waiting[i] = 1; if ( WIFEXITED(status) ) { continue; @@ -676,7 +717,7 @@ static void handle_zombie(struct sandbox *sb) STATUS("Last filename was: %s\n", sb->last_filename[i]); sb->n_processed++; - start_worker_process(sb, i); + start_worker_process(sb, i, 0, NULL); } } @@ -688,7 +729,7 @@ static void handle_zombie(struct sandbox *sb) void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, char *use_this_one_instead, - FILE *ofh) + FILE *ofh, int argc, char *argv[]) { int i; int allDone; @@ -703,17 +744,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } - sb->n_indexable = 0; sb->n_processed = 0; - sb->n_indexable_last_stats = 0; + sb->n_hadcrystals = 0; + sb->n_crystals = 0; sb->n_processed_last_stats = 0; + sb->n_hadcrystals_last_stats = 0; + sb->n_crystals_last_stats = 0; sb->t_last_stats = get_monotonic_seconds(); sb->n_proc = n_proc; - sb->ofh = ofh; sb->iargs = iargs; pthread_mutex_init(&sb->lock, NULL); + sb->ofh = ofh; sb->stream_pipe_read = calloc(n_proc, sizeof(int)); sb->stream_pipe_write = calloc(n_proc, sizeof(int)); if ( sb->stream_pipe_read == NULL ) { @@ -744,6 +787,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->result_fhs = calloc(n_proc, sizeof(FILE *)); sb->pids = calloc(n_proc, sizeof(pid_t)); sb->running = calloc(n_proc, sizeof(int)); + sb->waiting = calloc(n_proc, sizeof(int)); sb->fhs = calloc(sb->n_proc, sizeof(FILE *)); if ( sb->filename_pipes == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); @@ -761,6 +805,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Couldn't allocate memory for process flags.\n"); return; } + if ( sb->waiting == NULL ) { + ERROR("Couldn't allocate memory for process flags.\n"); + return; + } sb->last_filename = calloc(n_proc, sizeof(char *)); if ( sb->last_filename == NULL ) { @@ -796,7 +844,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, /* Fork the right number of times */ lock_sandbox(sb); for ( i=0; i<n_proc; i++ ) { - start_worker_process(sb, i); + start_worker_process(sb, i, argc, argv); } unlock_sandbox(sb); @@ -835,10 +883,9 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, r = select(fdmax+1, &fds, NULL, NULL, &tv); if ( r == -1 ) { - if ( errno != EINTR ) { - ERROR("select() failed: %s\n", strerror(errno)); - } - continue; + if ( errno == EINTR ) continue; + ERROR("select() failed: %s\n", strerror(errno)); + break; } if ( r == 0 ) continue; /* No progress this time. Try again */ @@ -886,8 +933,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Invalid result '%s'\n", results); } } else { - sb->n_indexable += atoi(results); + + int nc = atoi(results); + sb->n_crystals += nc; + if ( nc > 0 ) { + sb->n_hadcrystals++; + } sb->n_processed++; + } /* Send next filename */ @@ -920,14 +973,16 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, tNow = get_monotonic_seconds(); if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { - STATUS("%i out of %i indexed so far," - " %i out of %i since the last message.\n", - sb->n_indexable, sb->n_processed, - sb->n_indexable - sb->n_indexable_last_stats, + STATUS("%4i indexable out of %4i processed, " + "%4i crystals so far. " + "%4i images processed since the last message.\n", + sb->n_hadcrystals, sb->n_processed, + sb->n_crystals, sb->n_processed - sb->n_processed_last_stats); - sb->n_indexable_last_stats = sb->n_indexable; sb->n_processed_last_stats = sb->n_processed; + sb->n_hadcrystals_last_stats = sb->n_hadcrystals; + sb->n_crystals_last_stats = sb->n_crystals; sb->t_last_stats = tNow; } @@ -965,10 +1020,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->result_fhs); free(sb->pids); free(sb->running); + free(sb->waiting); - if ( ofh != stdout ) fclose(ofh); - - STATUS("There were %i images, of which %i could be indexed.\n", - sb->n_processed, sb->n_indexable); + STATUS("Final:" + " %i images processed, %i had crystals, %i crystals overall.\n", + sb->n_processed, sb->n_hadcrystals, sb->n_crystals); } diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 3851fb40..a454d867 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -31,6 +31,8 @@ * */ +#include "stream.h" + enum { PEAK_ZAEF, PEAK_HDF5, @@ -41,14 +43,12 @@ enum { struct index_args { UnitCell *cell; - int config_cmfilter; - int config_noisefilter; - int config_verbose; - int stream_flags; /* What goes into the output? */ - int config_satcorr; - int config_closer; - int config_insane; - int config_bgsub; + int cmfilter; + int noisefilter; + int verbose; + int satcorr; + int closer; + int bgsub; float threshold; float min_gradient; float min_snr; @@ -57,33 +57,25 @@ struct index_args IndexingMethod *indm; IndexingPrivate **ipriv; int peaks; /* Peak detection method */ - int cellr; float tols[4]; struct beam_params *beam; char *element; char *hdf5_peak_path; - double ir_inn; - double ir_mid; - double ir_out; + float ir_inn; + float ir_mid; + float ir_out; struct copy_hdf5_field *copyme; int integrate_saturated; int use_saturated; int no_revalidate; + int integrate_found; + int stream_peaks; + int stream_refls; + int res_cutoff; }; - - -/* Information about the indexing process for one pattern */ -struct pattern_args -{ - /* "Input" */ - char *filename; - - /* "Output" */ - int indexable; -}; - extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - char *use_this_one_instead, FILE *ofh); + char *use_this_one_instead, FILE *stream, + int argc, char *argv[]); diff --git a/src/indexamajig.c b/src/indexamajig.c index ae03986c..54ae0a5c 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -3,13 +3,13 @@ * * Index patterns, output hkl+intensity etc. * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -81,11 +81,8 @@ static void show_help(const char *s) " Default: indexamajig.stream\n" "\n" " --indexing=<methods> Use 'methods' for indexing. Provide one or more\n" -" methods separated by commas. Choose from:\n" -" none : no indexing (default)\n" -" dirax : invoke DirAx\n" -" mosflm : invoke MOSFLM (DPS)\n" -" reax : DPS algorithm with known unit cell\n" +" methods separated by commas.\n" +" See 'man indexamajig' for details.\n" " -g. --geometry=<file> Get detector geometry from file.\n" " -b, --beam=<file> Get beam parameters from file (provides nominal\n" " wavelength value if no per-shot value is found in\n" @@ -101,28 +98,8 @@ static void show_help(const char *s) " --hdf5-peaks=<p> Find peaks table in HDF5 file here.\n" " Default: /processing/hitfinder/peakinfo\n" "\n\n" -"You can control what information is included in the output stream using\n" -"' --record=<flag1>,<flag2>,<flag3>' and so on. Possible flags are:\n\n" -" integrated Include a list of reflection intensities, produced by\n" -" integrating around predicted peak locations.\n" -"\n" -" peaks Include peak locations and intensities from the peak\n" -" search.\n" -"\n" -" peaksifindexed As 'peaks', but only if the pattern could be indexed.\n" -"\n" -" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n" -"\n\n" -"The default is '--record=integrated'.\n" -"\n\n" "For more control over the process, you might need:\n\n" -" --cell-reduction=<m> Use <m> as the cell reduction method. Choose from:\n" -" none : no matching, just use the raw cell.\n" -" reduce : full cell reduction.\n" -" compare : match by at most changing the order of\n" -" the indices.\n" -" compare_ab : compare 'a' and 'b' lengths only.\n" -" --tolerance=<tol> Set the tolerances for cell reduction.\n" +" --tolerance=<tol> Set the tolerances for cell comparison.\n" " Default: 5,5,5,1.5.\n" " --filter-cm Perform common-mode noise subtraction on images\n" " before proceeding. Intensities will be extracted\n" @@ -134,7 +111,7 @@ static void show_help(const char *s) " --no-sat-corr Don't correct values of saturated peaks using a\n" " table included in the HDF5 file.\n" " --threshold=<n> Only accept peaks above <n> ADU. Default: 800.\n" -" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n" +" --min-gradient=<n> Minimum squared gradient for Zaefferer peak search.\n" " Default: 100,000.\n" " --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n" " Default: 5.\n" @@ -144,6 +121,8 @@ static void show_help(const char *s) "-e, --image=<element> Use this image from the HDF5 file.\n" " Example: /data/data0.\n" " Default: The first one found.\n" +" --res-cutoff Estimate the resolution limit of each pattern, and\n" +" don't integrate reflections further out.\n" "\n" "\nFor time-resolved stuff, you might want to use:\n\n" " --copy-hdf5-field <f> Copy the value of field <f> into the stream. You\n" @@ -157,8 +136,6 @@ static void show_help(const char *s) " --no-check-prefix Don't attempt to correct the --prefix.\n" " --closer-peak Don't integrate from the location of a nearby peak\n" " instead of the predicted spot. Don't use.\n" -" --insane Don't check that the reduced cell accounts for at\n" -" least 10%% of the located peaks.\n" " --no-bg-sub Don't subtract local background estimates from\n" " integrated intensities.\n" " --use-saturated During the initial peak search, don't reject\n" @@ -167,34 +144,15 @@ static void show_help(const char *s) " peaks which contain pixels above max_adu.\n" " --no-revalidate Don't re-integrate and check HDF5 peaks for\n" " validity.\n" +" --integrate-found Skip the spot prediction step, and just integrate\n" +" the intensities of the spots found by the initial\n" +" peak search.\n" +" --no-peaks-in-stream Do not record peak search results in the stream.\n" +" --no-refls-in-stream Do not record integrated reflections in the stream.\n" ); } -static int parse_cell_reduction(const char *scellr, int *err, - int *reduction_needs_cell) -{ - *err = 0; - if ( strcmp(scellr, "none") == 0 ) { - *reduction_needs_cell = 0; - return CELLR_NONE; - } else if ( strcmp(scellr, "reduce") == 0) { - *reduction_needs_cell = 1; - return CELLR_REDUCE; - } else if ( strcmp(scellr, "compare") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE; - } else if ( strcmp(scellr, "compare_ab") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE_AB; - } else { - *err = 1; - *reduction_needs_cell = 0; - return CELLR_NONE; - } -} - - int main(int argc, char *argv[]) { int c; @@ -203,87 +161,98 @@ int main(int argc, char *argv[]) FILE *fh; FILE *ofh; char *rval = NULL; - int config_noindex = 0; - int config_cmfilter = 0; - int config_noisefilter = 0; - int config_verbose = 0; - int config_satcorr = 1; int config_checkprefix = 1; - int config_closer = 0; - int config_insane = 0; - int config_bgsub = 1; int config_basename = 0; - float threshold = 800.0; - float min_gradient = 100000.0; - float min_snr = 5.0; - double min_int_snr = -INFINITY; - struct detector *det; - char *geometry = NULL; IndexingMethod *indm; IndexingPrivate **ipriv; - int indexer_needs_cell; - int reduction_needs_cell; char *indm_str = NULL; - UnitCell *cell; char *pdb = NULL; char *prefix = NULL; char *speaks = NULL; - char *scellr = NULL; char *toler = NULL; - float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */ - int cellr; - int peaks; int n_proc = 1; char *prepare_line; char prepare_filename[1024]; char *use_this_one_instead; struct index_args iargs; - struct beam_params *beam = NULL; - char *element = NULL; - double nominal_photon_energy; - int stream_flags = STREAM_INTEGRATED; - char *hdf5_peak_path = NULL; - struct copy_hdf5_field *copyme; char *intrad = NULL; - float ir_inn = 4.0; - float ir_mid = 5.0; - float ir_out = 7.0; - int integrate_saturated = 0; - int use_saturated = 0; - int no_revalidate = 0; - - copyme = new_copy_hdf5_field_list(); - if ( copyme == NULL ) { + + /* Defaults */ + iargs.cell = NULL; + iargs.cmfilter = 0; + iargs.noisefilter = 0; + iargs.verbose = 0; + iargs.satcorr = 1; + iargs.closer = 0; + iargs.bgsub = 1; + iargs.tols[0] = 5.0; + iargs.tols[1] = 5.0; + iargs.tols[2] = 5.0; + iargs.tols[3] = 1.5; + iargs.threshold = 800.0; + iargs.min_gradient = 100000.0; + iargs.min_snr = 5.0; + iargs.min_int_snr = -INFINITY; + iargs.det = NULL; + iargs.peaks = PEAK_ZAEF; + iargs.beam = NULL; + iargs.element = NULL; + iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); + iargs.copyme = NULL; + iargs.ir_inn = 4.0; + iargs.ir_mid = 5.0; + iargs.ir_out = 7.0; + iargs.use_saturated = 0; + iargs.integrate_saturated = 1; + iargs.no_revalidate = 0; + iargs.integrate_found = 0; + iargs.stream_peaks = 1; + iargs.stream_refls = 1; + iargs.res_cutoff = 0; + iargs.copyme = new_copy_hdf5_field_list(); + if ( iargs.copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); return 1; } + iargs.indm = NULL; /* No default */ + iargs.ipriv = NULL; /* No default */ /* Long options */ const struct option longopts[] = { + + /* Options with long and short versions */ {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, {"output", 1, NULL, 'o'}, - {"no-index", 0, &config_noindex, 1}, {"indexing", 1, NULL, 'z'}, {"geometry", 1, NULL, 'g'}, {"beam", 1, NULL, 'b'}, - {"filter-cm", 0, &config_cmfilter, 1}, - {"filter-noise", 0, &config_noisefilter, 1}, - {"verbose", 0, &config_verbose, 1}, {"pdb", 1, NULL, 'p'}, {"prefix", 1, NULL, 'x'}, - {"no-sat-corr", 0, &config_satcorr, 0}, - {"sat-corr", 0, &config_satcorr, 1}, /* Compat */ {"threshold", 1, NULL, 't'}, - {"no-check-prefix", 0, &config_checkprefix, 0}, - {"no-closer-peak", 0, &config_closer, 0}, - {"closer-peak", 0, &config_closer, 1}, - {"insane", 0, &config_insane, 1}, {"image", 1, NULL, 'e'}, - {"basename", 0, &config_basename, 1}, - {"bg-sub", 0, &config_bgsub, 1}, /* Compat */ - {"no-bg-sub", 0, &config_bgsub, 0}, + /* Long-only options with no arguments */ + {"filter-cm", 0, &iargs.cmfilter, 1}, + {"filter-noise", 0, &iargs.noisefilter, 1}, + {"verbose", 0, &iargs.verbose, 1}, + {"no-sat-corr", 0, &iargs.satcorr, 0}, + {"sat-corr", 0, &iargs.satcorr, 1}, + {"no-check-prefix", 0, &config_checkprefix, 0}, + {"no-closer-peak", 0, &iargs.closer, 0}, + {"closer-peak", 0, &iargs.closer, 1}, + {"basename", 0, &config_basename, 1}, + {"bg-sub", 0, &iargs.bgsub, 1}, + {"no-bg-sub", 0, &iargs.bgsub, 0}, + {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, + {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, + {"res-cutoff", 0, &iargs.res_cutoff, 1}, + {"integrate-saturated",0, &iargs.integrate_saturated,1}, + {"use-saturated", 0, &iargs.use_saturated, 1}, + {"no-revalidate", 0, &iargs.no_revalidate, 1}, + {"integrate-found", 0, &iargs.integrate_found, 1}, + + /* Long-only options with arguments */ {"peaks", 1, NULL, 2}, {"cell-reduction", 1, NULL, 3}, {"min-gradient", 1, NULL, 4}, @@ -297,9 +266,7 @@ int main(int argc, char *argv[]) {"min-integration-snr",1, NULL, 12}, {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, - {"integrate-saturated",0, &integrate_saturated,1}, - {"use-saturated",0, &use_saturated, 1}, - {"no-revalidate", 0, &no_revalidate, 1}, + {0, 0, NULL, 0} }; @@ -338,16 +305,21 @@ int main(int argc, char *argv[]) break; case 'g' : - geometry = strdup(optarg); + iargs.det = get_detector_geometry(optarg); + if ( iargs.det == NULL ) { + ERROR("Failed to read detector geometry from " + "'%s'\n", optarg); + return 1; + } break; case 't' : - threshold = strtof(optarg, NULL); + iargs.threshold = strtof(optarg, NULL); break; case 'b' : - beam = get_beam_parameters(optarg); - if ( beam == NULL ) { + iargs.beam = get_beam_parameters(optarg); + if ( iargs.beam == NULL ) { ERROR("Failed to load beam parameters" " from '%s'\n", optarg); return 1; @@ -355,7 +327,7 @@ int main(int argc, char *argv[]) break; case 'e' : - element = strdup(optarg); + iargs.element = strdup(optarg); break; case 2 : @@ -363,17 +335,24 @@ int main(int argc, char *argv[]) break; case 3 : - scellr = strdup(optarg); - break; + ERROR("The option '--cell-reduction' is no longer " + "used.\n" + "The complete indexing behaviour is now " + "controlled using '--indexing'.\n" + "See 'man indexamajig' for details of the " + "available methods.\n"); + return 1; case 4 : - min_gradient = strtof(optarg, NULL); + iargs.min_gradient = strtof(optarg, NULL); break; case 5 : - stream_flags = parse_stream_flags(optarg); - if ( stream_flags < 0 ) return 1; - break; + ERROR("The option '--record' is no longer used.\n" + "Use '--no-peaks-in-stream' and" + "'--no-refls-in-stream' if you need to control" + "the contents of the stream.\n"); + return 1; case 6 : case 7 : @@ -383,19 +362,20 @@ int main(int argc, char *argv[]) break; case 9 : - hdf5_peak_path = strdup(optarg); + free(iargs.hdf5_peak_path); + iargs.hdf5_peak_path = strdup(optarg); break; case 10 : - add_copy_hdf5_field(copyme, optarg); + add_copy_hdf5_field(iargs.copyme, optarg); break; case 11 : - min_snr = strtof(optarg, NULL); + iargs.min_snr = strtof(optarg, NULL); break; case 12 : - min_int_snr = strtof(optarg, NULL); + iargs.min_int_snr = strtof(optarg, NULL); break; case 13 : @@ -409,6 +389,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; @@ -431,30 +414,15 @@ int main(int argc, char *argv[]) } free(filename); - if ( outfile == NULL ) { - ofh = stdout; - } else { - ofh = fopen(outfile, "w"); - if ( ofh == NULL ) { - ERROR("Failed to open output file '%s'\n", outfile); - return 1; - } - free(outfile); - } - - if ( hdf5_peak_path == NULL ) { - hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); - } - if ( speaks == NULL ) { speaks = strdup("zaef"); STATUS("You didn't specify a peak detection method.\n"); STATUS("I'm using 'zaef' for you.\n"); } if ( strcmp(speaks, "zaef") == 0 ) { - peaks = PEAK_ZAEF; + iargs.peaks = PEAK_ZAEF; } else if ( strcmp(speaks, "hdf5") == 0 ) { - peaks = PEAK_HDF5; + iargs.peaks = PEAK_HDF5; } else { ERROR("Unrecognised peak detection method '%s'\n", speaks); return 1; @@ -474,57 +442,29 @@ int main(int argc, char *argv[]) return 1; } - if ( (indm_str == NULL) || - ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) { - STATUS("Not indexing anything.\n"); - indexer_needs_cell = 0; - reduction_needs_cell = 0; + if ( indm_str == NULL ) { + + STATUS("You didn't specify an indexing method, so I won't try " + " to index anything.\n" + "If that isn't what you wanted, re-run with" + " --indexing=<methods>.\n"); indm = NULL; - cellr = CELLR_NONE; + } else { - if ( indm_str == NULL ) { - STATUS("You didn't specify an indexing method, so I " - " won't try to index anything.\n" - "If that isn't what you wanted, re-run with" - " --indexing=<method>.\n"); - indm = NULL; - indexer_needs_cell = 0; - } else { - indm = build_indexer_list(indm_str, - &indexer_needs_cell); - if ( indm == NULL ) { - ERROR("Invalid indexer list '%s'\n", indm_str); - return 1; - } - free(indm_str); - } - reduction_needs_cell = 0; - if ( scellr == NULL ) { - STATUS("You didn't specify a cell reduction method, so" - " I'm going to use 'reduce'.\n"); - cellr = CELLR_REDUCE; - reduction_needs_cell = 1; - } else { - int err; - cellr = parse_cell_reduction(scellr, &err, - &reduction_needs_cell); - if ( err ) { - ERROR("Unrecognised cell reduction '%s'\n", - scellr); - return 1; - } - free(scellr); + indm = build_indexer_list(indm_str); + if ( indm == NULL ) { + ERROR("Invalid indexer list '%s'\n", indm_str); + return 1; } + free(indm_str); } - /* No indexing -> no reduction */ - if ( indm == NULL ) reduction_needs_cell = 0; - if ( toler != NULL ) { int ttt; ttt = sscanf(toler, "%f,%f,%f,%f", - &tols[0], &tols[1], &tols[2], &tols[3] ); + &iargs.tols[0], &iargs.tols[1], + &iargs.tols[2], &iargs.tols[3]); if ( ttt != 4 ) { ERROR("Invalid parameters for '--tolerance'\n"); return 1; @@ -534,7 +474,8 @@ int main(int argc, char *argv[]) if ( intrad != NULL ) { int r; - r = sscanf(intrad, "%f,%f,%f", &ir_inn, &ir_mid, &ir_out); + r = sscanf(intrad, "%f,%f,%f", + &iargs.ir_inn, &iargs.ir_mid, &iargs.ir_out); if ( r != 3 ) { ERROR("Invalid parameters for '--int-radius'\n"); return 1; @@ -546,43 +487,38 @@ int main(int argc, char *argv[]) " probably not appropriate for your patterns.\n"); } - if ( geometry == NULL ) { - ERROR("You need to specify a geometry file with --geometry\n"); + if ( iargs.det == NULL ) { + ERROR("You need to provide a geometry file (please read the" + " manual for more details).\n"); return 1; } - det = get_detector_geometry(geometry); - if ( det == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", geometry); + if ( iargs.beam == NULL ) { + ERROR("You need to provide a beam parameters file (please read" + " the manual for more details).\n"); return 1; } - free(geometry); if ( pdb != NULL ) { - cell = load_cell_from_pdb(pdb); - if ( cell == NULL ) { + iargs.cell = load_cell_from_pdb(pdb); + if ( iargs.cell == NULL ) { ERROR("Couldn't read unit cell (from %s)\n", pdb); return 1; } free(pdb); - cell_print(cell); + STATUS("This is what I understood your unit cell to be:\n"); + cell_print(iargs.cell); } else { STATUS("No unit cell given.\n"); - cell = NULL; + iargs.cell = NULL; } - write_stream_header(ofh, argc, argv); - - if ( beam != NULL ) { - nominal_photon_energy = beam->photon_energy; - } else { - STATUS("No beam parameters file was given, so I'm taking the" - " nominal photon energy to be 2 keV.\n"); - ERROR("I'm also going to assume 1 ADU per photon, which is"); - ERROR(" almost certainly wrong. Peak sigmas will be" - " incorrect.\n"); - nominal_photon_energy = 2000.0; + ofh = fopen(outfile, "w"); + if ( ofh == NULL ) { + ERROR("Failed to open stream '%s'\n", outfile); + return 1; } + free(outfile); /* Get first filename and use it to set up the indexing */ prepare_line = malloc(1024); @@ -604,8 +540,8 @@ int main(int argc, char *argv[]) /* Prepare the indexer */ if ( indm != NULL ) { - ipriv = prepare_indexing(indm, cell, prepare_filename, det, - nominal_photon_energy); + ipriv = prepare_indexing(indm, iargs.cell, prepare_filename, + iargs.det, iargs.beam, iargs.tols); if ( ipriv == NULL ) { ERROR("Failed to prepare indexing.\n"); return 1; @@ -616,42 +552,11 @@ int main(int argc, char *argv[]) gsl_set_error_handler_off(); - /* Static worker args */ - iargs.cell = cell; - iargs.config_cmfilter = config_cmfilter; - iargs.config_noisefilter = config_noisefilter; - iargs.config_verbose = config_verbose; - iargs.config_satcorr = config_satcorr; - iargs.config_closer = config_closer; - iargs.config_insane = config_insane; - iargs.config_bgsub = config_bgsub; - iargs.cellr = cellr; - iargs.tols[0] = tols[0]; - iargs.tols[1] = tols[1]; - iargs.tols[2] = tols[2]; - iargs.tols[3] = tols[3]; - iargs.threshold = threshold; - iargs.min_gradient = min_gradient; - iargs.min_snr = min_snr; - iargs.min_int_snr = min_int_snr; - iargs.det = det; iargs.indm = indm; iargs.ipriv = ipriv; - iargs.peaks = peaks; - iargs.beam = beam; - iargs.element = element; - iargs.stream_flags = stream_flags; - iargs.hdf5_peak_path = hdf5_peak_path; - iargs.copyme = copyme; - iargs.ir_inn = ir_inn; - iargs.ir_mid = ir_mid; - iargs.ir_out = ir_out; - iargs.use_saturated = use_saturated; - iargs.integrate_saturated = integrate_saturated; - iargs.no_revalidate = no_revalidate; create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - use_this_one_instead, ofh); + use_this_one_instead, ofh, argc, argv); free(prefix); diff --git a/src/partial_sim.c b/src/partial_sim.c index e9a7003e..e48c51fa 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -54,11 +54,12 @@ #define NBINS 50 -static void mess_up_cell(UnitCell *cell, double cnoise) +static void mess_up_cell(Crystal *cr, double cnoise) { double ax, ay, az; double bx, by, bz; double cx, cy, cz; + UnitCell *cell = crystal_get_cell(cr); //STATUS("Real:\n"); //cell_print(cell); @@ -82,17 +83,18 @@ static void mess_up_cell(UnitCell *cell, double cnoise) /* For each reflection in "partial", fill in what the intensity would be * according to "full" */ -static void calculate_partials(RefList *partial, double osf, +static void calculate_partials(Crystal *cr, RefList *full, const SymOpList *sym, int random_intensities, pthread_mutex_t *full_lock, unsigned long int *n_ref, double *p_hist, - double *p_max, double max_q, UnitCell *cell) + double *p_max, double max_q) { Reflection *refl; RefListIterator *iter; + double res; - for ( refl = first_refl(partial, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -137,16 +139,17 @@ static void calculate_partials(RefList *partial, double osf, } } - Ip = osf * p * If; + Ip = crystal_get_osf(cr) * p * If; - bin = NBINS*2.0*resolution(cell, h, k, l) / max_q; + res = resolution(crystal_get_cell(cr), h, k, l); + bin = NBINS*2.0*res/max_q; if ( (bin < NBINS) && (bin>=0) ) { p_hist[bin] += p; n_ref[bin]++; if ( p > p_max[bin] ) p_max[bin] = p; } else { STATUS("Reflection out of histogram range: %e %i %f\n", - resolution(cell, h, k, l), bin, p); + res, bin, p); } Ip = gaussian_noise(Ip, 100.0); @@ -206,13 +209,14 @@ struct queue_args unsigned long int n_ref[NBINS]; double p_max[NBINS]; - FILE *stream; + Stream *stream; }; struct worker_args { struct queue_args *qargs; + Crystal *crystal; struct image image; /* Histogram for this image */ @@ -243,21 +247,31 @@ static void *create_job(void *vqargs) static void run_job(void *vwargs, int cookie) { - double osf; struct quaternion orientation; struct worker_args *wargs = vwargs; struct queue_args *qargs = wargs->qargs; int i; + Crystal *cr; + RefList *reflections; - osf = gaussian_noise(1.0, 0.3); + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to create crystal.\n"); + return; + } + wargs->crystal = cr; + crystal_set_image(cr, &wargs->image); + + crystal_set_osf(cr, gaussian_noise(1.0, 0.3)); + crystal_set_profile_radius(cr, wargs->image.beam->profile_radius); /* Set up a random orientation */ orientation = random_quaternion(); - wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation); + crystal_set_cell(cr, cell_rotate(qargs->cell, orientation)); snprintf(wargs->image.filename, 255, "dummy.h5"); - wargs->image.reflections = find_intersections(&wargs->image, - wargs->image.indexed_cell); + reflections = find_intersections(&wargs->image, cr); + crystal_set_reflections(cr, reflections); for ( i=0; i<NBINS; i++ ) { wargs->n_ref[i] = 0; @@ -265,14 +279,16 @@ static void run_job(void *vwargs, int cookie) wargs->p_max[i] = 0.0; } - calculate_partials(wargs->image.reflections, osf, qargs->full, + calculate_partials(cr, qargs->full, qargs->sym, qargs->random_intensities, &qargs->full_lock, wargs->n_ref, wargs->p_hist, wargs->p_max, - qargs->max_q, wargs->image.indexed_cell); + qargs->max_q); /* Give a slightly incorrect cell in the stream */ - mess_up_cell(wargs->image.indexed_cell, qargs->cnoise); + mess_up_cell(cr, qargs->cnoise); + + image_add_crystal(&wargs->image, cr); } @@ -282,7 +298,7 @@ static void finalise_job(void *vqargs, void *vwargs) struct queue_args *qargs = vqargs; int i; - write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED); + write_chunk(qargs->stream, &wargs->image, NULL, 0, 1); for ( i=0; i<NBINS; i++ ) { qargs->n_ref[i] += wargs->n_ref[i]; @@ -295,8 +311,7 @@ static void finalise_job(void *vqargs, void *vwargs) qargs->n_done++; progress_bar(qargs->n_done, qargs->n_to_do, "Simulating"); - reflist_free(wargs->image.reflections); - cell_free(wargs->image.indexed_cell); + crystal_free(wargs->crystal); free(wargs); } @@ -315,7 +330,7 @@ int main(int argc, char *argv[]) char *sym_str = NULL; SymOpList *sym; UnitCell *cell = NULL; - FILE *ofh; + Stream *stream; int n = 2; int random_intensities = 0; char *save_file = NULL; @@ -404,6 +419,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; @@ -493,13 +511,12 @@ int main(int argc, char *argv[]) ERROR("You must give a filename for the output.\n"); return 1; } - ofh = fopen(output_file, "w"); - if ( ofh == NULL ) { + stream = open_stream_for_write(output_file); + if ( stream == NULL ) { ERROR("Couldn't open output file '%s'\n", output_file); return 1; } free(output_file); - write_stream_header(ofh, argc, argv); image.det = det; image.width = det->max_fs; @@ -508,9 +525,11 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy)); image.div = beam->divergence; image.bw = beam->bandwidth; - image.profile_radius = beam->profile_radius; + image.beam = beam; image.filename = malloc(256); image.copyme = NULL; + image.crystals = NULL; + image.n_crystals = 0; if ( random_intensities ) { full = reflist_new(); @@ -525,7 +544,7 @@ int main(int argc, char *argv[]) qargs.random_intensities = random_intensities; qargs.cell = cell; qargs.template_image = ℑ - qargs.stream = ofh; + qargs.stream = stream; qargs.cnoise = cnoise; qargs.max_q = largest_q(&image); @@ -571,7 +590,7 @@ int main(int argc, char *argv[]) } - fclose(ofh); + close_stream(stream); cell_free(cell); free_detector_geometry(det); free(beam); diff --git a/src/partialator.c b/src/partialator.c index 5d373f83..8b51f826 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -87,16 +87,16 @@ static void show_help(const char *s) struct refine_args { RefList *full; - struct image *image; + Crystal *crystal; }; struct queue_args { - int n; + int n_started; int n_done; - int n_total_patterns; - struct image *images; + Crystal **crystals; + int n_crystals; struct refine_args task_defaults; }; @@ -104,10 +104,9 @@ struct queue_args static void refine_image(void *task, int id) { struct refine_args *pargs = task; - struct image *image = pargs->image; - image->id = id; + Crystal *cr = pargs->crystal; - pr_refine(image, pargs->full); + pr_refine(cr, pargs->full); } @@ -119,9 +118,9 @@ static void *get_image(void *vqargs) task = malloc(sizeof(struct refine_args)); memcpy(task, &qargs->task_defaults, sizeof(struct refine_args)); - task->image = &qargs->images[qargs->n]; + task->crystal = qargs->crystals[qargs->n_started]; - qargs->n++; + qargs->n_started++; return task; } @@ -133,12 +132,12 @@ static void done_image(void *vqargs, void *task) qargs->n_done++; - progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining"); + progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); free(task); } -static void refine_all(struct image *images, int n_total_patterns, +static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, RefList *full, int nthreads) { @@ -146,19 +145,19 @@ static void refine_all(struct image *images, int n_total_patterns, struct queue_args qargs; task_defaults.full = full; - task_defaults.image = NULL; + task_defaults.crystal = NULL; qargs.task_defaults = task_defaults; - qargs.n = 0; + qargs.n_started = 0; qargs.n_done = 0; - qargs.n_total_patterns = n_total_patterns; - qargs.images = images; + qargs.n_crystals = n_crystals; + qargs.crystals = crystals; /* Don't have threads which are doing nothing */ - if ( n_total_patterns < nthreads ) nthreads = n_total_patterns; + if ( n_crystals < nthreads ) nthreads = n_crystals; run_threads(nthreads, refine_image, get_image, done_image, - &qargs, n_total_patterns, 0, 0, 0); + &qargs, n_crystals, 0, 0, 0); } @@ -201,13 +200,14 @@ static int select_scalable_reflections(RefList *list, RefList *reference) } -static void select_reflections_for_refinement(struct image *images, int n, +static void select_reflections_for_refinement(Crystal **crystals, int n, RefList *full, int have_reference) { int i; for ( i=0; i<n; i++ ) { + RefList *reflist; Reflection *refl; RefListIterator *iter; int n_acc = 0; @@ -216,9 +216,10 @@ static void select_reflections_for_refinement(struct image *images, int n, int n_fewmatch = 0; int n_ref = 0; - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(crystals[i]) ) continue; - for ( refl = first_refl(images[i].reflections, &iter); + reflist = crystal_get_reflections(crystals[i]); + for ( refl = first_refl(reflist, &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -272,6 +273,20 @@ static void select_reflections_for_refinement(struct image *images, int n, } +static void display_progress(int n_images, int n_crystals) +{ + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%i images loaded, %i crystals.", + n_images, n_crystals); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + int main(int argc, char *argv[]) { int c; @@ -280,16 +295,15 @@ int main(int argc, char *argv[]) char *geomfile = NULL; char *sym_str = NULL; SymOpList *sym; - FILE *fh; int nthreads = 1; struct detector *det; int i; - int n_total_patterns; struct image *images; int n_iter = 10; struct beam_params *beam = NULL; RefList *full; - int n_usable_patterns = 0; + int n_images = 0; + int n_crystals = 0; int nobs; char *reference_file = NULL; RefList *reference = NULL; @@ -298,6 +312,8 @@ int main(int argc, char *argv[]) char cmdline[1024]; SRContext *sr; int noscale = 0; + Stream *st; + Crystal **crystals; /* Long options */ const struct option longopts[] = { @@ -370,6 +386,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; @@ -383,19 +402,15 @@ int main(int argc, char *argv[]) return 1; } - /* Sanitise input filename and open */ if ( infile == NULL ) { infile = strdup("-"); } - if ( strcmp(infile, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(infile, "r"); - } - if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", infile); + st = open_stream_for_read(infile); + if ( st == NULL ) { + ERROR("Failed to open input stream '%s'\n", infile); return 1; } + free(infile); /* Sanitise output filename */ if ( outfile == NULL ) { @@ -435,86 +450,96 @@ int main(int argc, char *argv[]) } - n_total_patterns = count_patterns(fh); - if ( n_total_patterns == 0 ) { - ERROR("No patterns to process.\n"); - return 1; - } - STATUS("There are %i patterns to process\n", n_total_patterns); - gsl_set_error_handler_off(); - images = malloc(n_total_patterns * sizeof(struct image)); - if ( images == NULL ) { - ERROR("Couldn't allocate memory for images.\n"); - return 1; - } - /* Fill in what we know about the images so far */ - rewind(fh); nobs = 0; - for ( i=0; i<n_total_patterns; i++ ) { + n_images = 0; + n_crystals = 0; + images = NULL; + crystals = NULL; + do { RefList *as; + int i; + struct image *images_new; struct image *cur; - cur = &images[n_usable_patterns]; + images_new = realloc(images, (n_images+1)*sizeof(struct image)); + if ( images_new == NULL ) { + ERROR("Failed to allocate memory for image list.\n"); + return 1; + } + images = images_new; + cur = &images[n_images]; cur->det = det; - - if ( read_chunk(fh, cur) != 0 ) { - /* Should not happen, because we counted the patterns - * earlier. */ - ERROR("Failed to read chunk from the input stream.\n"); - return 1; + if ( read_chunk(st, cur) != 0 ) { + break; } /* Won't be needing this, if it exists */ image_feature_list_free(cur->features); cur->features = NULL; - - /* "n_usable_patterns" will not be incremented in this case */ - if ( cur->indexed_cell == NULL ) continue; - - /* Fill in initial estimates of stuff */ cur->div = beam->divergence; cur->bw = beam->bandwidth; cur->width = det->max_fs; cur->height = det->max_ss; - cur->osf = 1.0; - cur->profile_radius = beam->profile_radius; - cur->pr_dud = 0; - - /* Muppet proofing */ cur->data = NULL; cur->flags = NULL; cur->beam = NULL; - /* This is the raw list of reflections */ - as = asymmetric_indices(cur->reflections, sym); - reflist_free(cur->reflections); - cur->reflections = as; + n_images++; + + for ( i=0; i<cur->n_crystals; i++ ) { + + Crystal *cr; + Crystal **crystals_new; - update_partialities(cur); + crystals_new = realloc(crystals, + (n_crystals+1)*sizeof(Crystal *)); + if ( crystals_new == NULL ) { + ERROR("Failed to allocate memory for crystal " + "list.\n"); + return 1; + } + crystals = crystals_new; + crystals[n_crystals] = cur->crystals[i]; + cr = crystals[n_crystals]; + crystal_set_image(cr, cur); - nobs += select_scalable_reflections(cur->reflections, - reference); + /* Fill in initial estimates of stuff */ + crystal_set_osf(cr, 1.0); + crystal_set_profile_radius(cr, beam->profile_radius); + crystal_set_user_flag(cr, 0); - progress_bar(i, n_total_patterns-1, "Loading pattern data"); - n_usable_patterns++; + /* This is the raw list of reflections */ + as = asymmetric_indices(crystal_get_reflections(cr), + sym); + crystal_set_reflections(cr, as); + update_partialities(cr); - } - fclose(fh); + nobs += select_scalable_reflections(as, reference); + + n_crystals++; + + } + + display_progress(n_images, n_crystals); + + } while ( 1 ); + + close_stream(st); /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); + STATUS("\nPerforming initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); - full = scale_intensities(images, n_usable_patterns, reference, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf", + sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); - sr_iteration(sr, 0, images, n_usable_patterns, full); + sr_iteration(sr, 0, crystals, n_crystals, full); /* Iterate */ for ( i=0; i<n_iter; i++ ) { @@ -531,34 +556,34 @@ int main(int argc, char *argv[]) } /* Refine the geometry of all patterns to get the best fit */ - select_reflections_for_refinement(images, n_usable_patterns, + select_reflections_for_refinement(crystals, n_crystals, comp, have_reference); - refine_all(images, n_usable_patterns, det, comp, nthreads); + refine_all(crystals, n_crystals, det, comp, nthreads); nobs = 0; - for ( j=0; j<n_usable_patterns; j++ ) { + for ( j=0; j<n_crystals; j++ ) { - struct image *cur = &images[j]; + Crystal *cr = crystals[j]; + RefList *rf = crystal_get_reflections(cr); - nobs += select_scalable_reflections(cur->reflections, - reference); + nobs += select_scalable_reflections(rf, reference); } /* Re-estimate all the full intensities */ reflist_free(full); - full = scale_intensities(images, n_usable_patterns, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr_iteration(sr, i+1, images, n_usable_patterns, full); + sr_iteration(sr, i+1, crystals, n_crystals, full); } sr_finish(sr); n_dud = 0; - for ( i=0; i<n_usable_patterns; i++ ) { - if ( images[i].pr_dud ) n_dud++; + for ( i=0; i<n_crystals; i++ ) { + if ( crystal_get_user_flag(crystals[i]) ) n_dud++; } STATUS("%i images could not be refined on the last cycle.\n", n_dud); @@ -566,19 +591,19 @@ int main(int argc, char *argv[]) write_reflist(outfile, full); /* Clean up */ - for ( i=0; i<n_usable_patterns; i++ ) { - reflist_free(images[i].reflections); + for ( i=0; i<n_crystals; i++ ) { + crystal_free(crystals[i]); } reflist_free(full); free(sym); free(outfile); free_detector_geometry(det); free(beam); + free(crystals); if ( reference != NULL ) { reflist_free(reference); } - for ( i=0; i<n_usable_patterns; i++ ) { - cell_free(images[i].indexed_cell); + for ( i=0; i<n_images; i++ ) { free(images[i].filename); } free(images); diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 89413a07..606a173c 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -369,6 +369,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; diff --git a/src/post-refinement.c b/src/post-refinement.c index 7410d931..1439b148 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, Reflection *refl, double r) +double gradient(Crystal *cr, int k, Reflection *refl) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -103,17 +103,19 @@ double gradient(struct image *image, int k, Reflection *refl, double r) int clamp_low, clamp_high; double klow, khigh; double gr; + struct image *image = crystal_get_image(cr); + double r = crystal_get_profile_radius(cr); get_symmetric_indices(refl, &hs, &ks, &ls); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); xl = hs*asx + ks*bsx + ls*csx; yl = hs*asy + ks*bsy + ls*csy; zl = hs*asz + ks*bsz + ls*csz; - ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); @@ -237,8 +239,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) /* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(struct image *image, int k, double shift) +static void apply_shift(Crystal *cr, int k, double shift) { + double t; + struct image *image = crystal_get_image(cr); + switch ( k ) { case REF_DIV : @@ -251,7 +256,9 @@ static void apply_shift(struct image *image, int k, double shift) break; case REF_R : - image->profile_radius += shift; + t = crystal_get_profile_radius(cr); + t += shift; + crystal_set_profile_radius(cr, t); break; case REF_ASX : @@ -263,7 +270,7 @@ static void apply_shift(struct image *image, int k, double shift) case REF_CSX : case REF_CSY : case REF_CSZ : - apply_cell_shift(image->indexed_cell, k, shift); + apply_cell_shift(crystal_get_cell(cr), k, shift); break; default : @@ -357,7 +364,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(struct image *image, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full) { gsl_matrix *M; gsl_vector *v; @@ -369,7 +376,7 @@ static double pr_iterate(struct image *image, const RefList *full) double max_shift; int nref = 0; - reflections = image->reflections; + reflections = crystal_get_reflections(cr); M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); @@ -387,6 +394,7 @@ static double pr_iterate(struct image *image, const RefList *full) double p; Reflection *match; double gradients[NUM_PARAMS]; + const double osf = crystal_get_osf(cr); if ( !get_refinable(refl) ) continue; @@ -398,7 +406,7 @@ static double pr_iterate(struct image *image, const RefList *full) " is it marked as refinable?\n", ha, ka, la); continue; } - I_full = image->osf * get_intensity(match); + I_full = osf * get_intensity(match); /* Actual measurement of this reflection from this pattern? */ I_partial = get_intensity(refl); @@ -412,7 +420,7 @@ static double pr_iterate(struct image *image, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; k<NUM_PARAMS; k++ ) { double gr; - gr = gradient(image, k, refl, image->profile_radius); + gr = gradient(cr, k, refl); gradients[k] = gr; } @@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full) if ( g > k ) continue; M_c = gradients[g] * gradients[k]; - M_c *= w * pow(image->osf * I_full, 2.0); + M_c *= w * pow(osf * I_full, 2.0); M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full) //STATUS("%i reflections went into the equations.\n", nref); if ( nref == 0 ) { - image->pr_dud = 1; + crystal_set_user_flag(cr, 1); gsl_matrix_free(M); gsl_vector_free(v); return 0.0; @@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full) for ( param=0; param<NUM_PARAMS; param++ ) { double shift = gsl_vector_get(shifts, param); - apply_shift(image, param, shift); + apply_shift(cr, param, shift); //STATUS("Shift %i: %e\n", param, shift); if ( fabs(shift) > max_shift ) max_shift = fabs(shift); } @@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full) } -static double guide_dev(struct image *image, const RefList *full) +static double guide_dev(Crystal *cr, const RefList *full) { double dev = 0.0; @@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full) Reflection *refl; RefListIterator *iter; - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full) * scale_intensities() might not yet have been called, so the * full version may not have been calculated yet. */ - G = image->osf; + G = crystal_get_osf(cr); p = get_partiality(refl); I_partial = get_intensity(refl); I_full = get_intensity(full_version); @@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full) } -void pr_refine(struct image *image, const RefList *full) +void pr_refine(Crystal *cr, const RefList *full) { double max_shift, dev; int i; const int verbose = 0; if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("\n"); /* Deal with progress bar */ STATUS("Before iteration: dev = %10.5e\n", dev); } i = 0; - image->pr_dud = 0; + crystal_set_user_flag(cr, 0); do { double asx, asy, asz; @@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full) double csx, csy, csz; double dev; - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(image, full); + max_shift = pr_iterate(cr, full); - update_partialities(image); + update_partialities(cr); if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("PR Iteration %2i: max shift = %10.2f" " dev = %10.5e\n", i+1, max_shift, dev); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 2223dcdf..fe171882 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -39,6 +39,7 @@ #include "image.h" #include "utils.h" +#include "crystal.h" /* Refineable parameters */ @@ -58,10 +59,10 @@ enum { }; -extern void pr_refine(struct image *image, const RefList *full); +extern void pr_refine(Crystal *cr, const RefList *full); /* Exported so it can be poked by tests/pr_gradient_check */ -extern double gradient(struct image *image, int k, Reflection *refl, double r); +extern double gradient(Crystal *cr, int k, Reflection *refl); #endif /* POST_REFINEMENT_H */ diff --git a/src/process_hkl.c b/src/process_hkl.c index cf690e9b..a8abc2e8 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -3,12 +3,12 @@ * * Assemble and process FEL Bragg intensities * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White <taw@physics.org> + * 2009-2013 Thomas White <taw@physics.org> * 2011 Andrew Martin <andrew.martin@desy.de> * 2012 Lorenzo Galli <lorenzo.galli@desy.de> * @@ -48,6 +48,8 @@ #include "stream.h" #include "reflist.h" #include "image.h" +#include "crystal.h" +#include "thread-pool.h" static void show_help(const char *s) @@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new, } -static int merge_pattern(RefList *model, struct image *new, RefList *reference, - const SymOpList *sym, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, int *hist_n, - int config_nopolar) +static int merge_crystal(RefList *model, struct image *image, Crystal *cr, + RefList *reference, const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, int *hist_n, + int config_nopolar) { Reflection *refl; RefListIterator *iter; @@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, double csx, csy, csz; if ( reference != NULL ) { - scale = scale_intensities(reference, new->reflections, sym); + scale = scale_intensities(reference, + crystal_get_reflections(cr), sym); } else { scale = 1.0; } if ( isnan(scale) ) return 1; - cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); - for ( refl = first_refl(new->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - ool = 1.0 / new->lambda; + ool = 1.0 / image->lambda; tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); phi = atan2(yl, xl); pa = pow(sin(phi)*sin(tt), 2.0); @@ -260,60 +263,67 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, } -static void merge_all(FILE *fh, RefList *model, RefList *reference, - int config_startafter, int config_stopafter, - const SymOpList *sym, - int n_total_patterns, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, - int *hist_i, int config_nopolar, int min_measurements) +static void display_progress(int n_images, int n_crystals, int n_crystals_used) +{ + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used.", + n_images, n_crystals, n_crystals_used); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + + +static int merge_all(Stream *st, RefList *model, RefList *reference, + const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, + int *hist_i, int config_nopolar, int min_measurements) { int rval; - int n_patterns = 0; - int n_used = 0; + int n_images = 0; + int n_crystals = 0; + int n_crystals_used = 0; Reflection *refl; RefListIterator *iter; - if ( skip_some_files(fh, config_startafter) ) { - ERROR("Failed to skip first %i files.\n", config_startafter); - return; - } - do { struct image image; + int i; image.det = NULL; /* Get data from next chunk */ - rval = read_chunk(fh, &image); + rval = read_chunk(st, &image); if ( rval ) break; - n_patterns++; + n_images++; - if ( (image.reflections != NULL) && (image.indexed_cell) ) { + for ( i=0; i<image.n_crystals; i++ ) { int r; + Crystal *cr = image.crystals[i]; - r = merge_pattern(model, &image, reference, sym, + n_crystals++; + r = merge_crystal(model, &image, cr, reference, sym, hist_vals, hist_h, hist_k, hist_l, hist_i, config_nopolar); - if ( r == 0 ) n_used++; + if ( r == 0 ) n_crystals_used++; + + crystal_free(cr); } free(image.filename); - reflist_free(image.reflections); image_feature_list_free(image.features); - cell_free(image.indexed_cell); - progress_bar(n_patterns, n_total_patterns-config_startafter, - "Merging"); - - if ( config_stopafter ) { - if ( n_patterns == config_stopafter ) break; - } + display_progress(n_images, n_crystals, n_crystals_used); } while ( rval == 0 ); @@ -339,7 +349,7 @@ static void merge_all(FILE *fh, RefList *model, RefList *reference, } - STATUS("%i of the patterns could be used.\n", n_used); + return 0; } @@ -348,14 +358,11 @@ int main(int argc, char *argv[]) int c; char *filename = NULL; char *output = NULL; - FILE *fh; + Stream *st; RefList *model; int config_maxonly = 0; - int config_startafter = 0; - int config_stopafter = 0; int config_sum = 0; int config_scale = 0; - unsigned int n_total_patterns; char *sym_str = NULL; SymOpList *sym; char *histo = NULL; @@ -369,6 +376,7 @@ int main(int argc, char *argv[]) int config_nopolar = 0; char *rval; int min_measurements = 2; + int r; /* Long options */ const struct option longopts[] = { @@ -409,11 +417,11 @@ int main(int argc, char *argv[]) break; case 's' : - config_stopafter = atoi(optarg); + ERROR("The option '--stop-after' no longer works.\n"); break; case 'f' : - config_startafter = atoi(optarg); + ERROR("The option '--start-after' no longer works.\n"); break; case 'y' : @@ -437,6 +445,9 @@ int main(int argc, char *argv[]) } break; + case '?' : + break; + case 0 : break; @@ -462,25 +473,11 @@ int main(int argc, char *argv[]) free(sym_str); /* Open the data stream */ - if ( strcmp(filename, "-") == 0 ) { - fh = stdin; - } else { - fh = fopen(filename, "r"); - } - free(filename); - if ( fh == NULL ) { - ERROR("Failed to open input file\n"); - return 1; - } - - /* Count the number of patterns in the file */ - n_total_patterns = count_patterns(fh); - if ( n_total_patterns == 0 ) { - ERROR("No patterns to process.\n"); + st = open_stream_for_read(filename); + if ( st == NULL ) { + ERROR("Failed to open stream.\n"); return 1; } - STATUS("There are %i patterns to process\n", n_total_patterns); - rewind(fh); model = reflist_new(); @@ -494,7 +491,8 @@ int main(int argc, char *argv[]) return 1; } - space_for_hist = n_total_patterns * num_equivs(sym, NULL); + /* FIXME: This array must grow as necessary. */ + space_for_hist = 0 * num_equivs(sym, NULL); hist_vals = malloc(space_for_hist * sizeof(double)); free(histo); STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l); @@ -526,40 +524,46 @@ int main(int argc, char *argv[]) } hist_i = 0; - merge_all(fh, model, NULL, config_startafter, config_stopafter, - sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l, - &hist_i, config_nopolar, min_measurements); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); + r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l, + &hist_i, config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); return 1; } - rewind(fh); if ( config_scale ) { RefList *reference; - STATUS("Extra pass for scaling...\n"); + if ( rewind_stream(st) ) { - reference = copy_reflist(model); + ERROR("Couldn't rewind stream - scaling cannot be " + "performed.\n"); - reflist_free(model); - model = reflist_new(); + } else { - rewind(fh); + int r; - merge_all(fh, model, reference, - config_startafter, config_stopafter, sym, - n_total_patterns, - hist_vals, hist_h, hist_k, hist_l, &hist_i, - config_nopolar, min_measurements); + STATUS("Extra pass for scaling...\n"); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); - return 1; - } + reference = copy_reflist(model); + + reflist_free(model); + model = reflist_new(); + + r = merge_all(st, model, reference, sym, + hist_vals, hist_h, hist_k, hist_l, &hist_i, + config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); + return 1; + } + + reflist_free(reference); - reflist_free(reference); + } } @@ -576,7 +580,7 @@ int main(int argc, char *argv[]) write_reflist(output, model); - fclose(fh); + close_stream(st); free(sym); reflist_free(model); diff --git a/src/render_hkl.c b/src/render_hkl.c index dbcdba43..142438e7 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -704,6 +704,9 @@ int main(int argc, char *argv[]) case 0 : break; + case '?' : + break; + default : ERROR("Unhandled option '%c'\n", c); break; diff --git a/src/scaling-report.c b/src/scaling-report.c index ba425b96..b161a01e 100644 --- a/src/scaling-report.c +++ b/src/scaling-report.c @@ -3,11 +3,11 @@ * * Write a nice PDF of scaling parameters * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height, } -static void partiality_graph(cairo_t *cr, const struct image *images, int n, +static void partiality_graph(cairo_t *cr, Crystal **crystals, int n, RefList *full) { const double g_width = 200.0; @@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, num_nondud = 0; for ( i=0; i<n; i++ ) { - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(crystals[i]) ) continue; num_nondud++; } @@ -229,10 +229,13 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, Reflection *refl; RefListIterator *iter; + Crystal *cryst; - if ( images[i].pr_dud ) continue; + cryst = crystals[i]; - for ( refl = first_refl(images[i].reflections, &iter); + if ( crystal_get_user_flag(cryst) ) continue; + + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -249,7 +252,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, if ( get_redundancy(f) < 2 ) continue; Ipart = get_intensity(refl); - Ifull = images[i].osf * get_intensity(f); + Ifull = crystal_get_osf(cryst) * get_intensity(f); //if ( Ifull < 10 ) continue; /* FIXME: Ugh */ @@ -267,7 +270,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, double esd_pobs, esd_Ip, esd_If; esd_Ip = get_esd_intensity(refl); esd_If = get_esd_intensity(f); - esd_If *= images[i].osf; + esd_If *= crystal_get_osf(cryst); esd_pobs = pow(esd_Ip/Ipart, 2.0); esd_pobs += pow(esd_If/Ifull, 2.0); esd_pobs = sqrt(esd_pobs); @@ -313,8 +316,8 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, } -static void partiality_histogram(cairo_t *cr, const struct image *images, - int n, RefList *full, int calc, int backwards) +static void partiality_histogram(cairo_t *cr, Crystal **crystals, int n, + RefList *full, int calc, int backwards) { int f_max; int i, b; @@ -344,10 +347,11 @@ static void partiality_histogram(cairo_t *cr, const struct image *images, Reflection *refl; RefListIterator *iter; + Crystal *cryst = crystals[i]; - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(cryst) ) continue; - for ( refl = first_refl(images[i].reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -364,7 +368,7 @@ static void partiality_histogram(cairo_t *cr, const struct image *images, Ipart = get_intensity(refl); Ifull = get_intensity(f); - pobs = Ipart/(images[i].osf*Ifull); + pobs = Ipart/(crystal_get_osf(cryst)*Ifull); pcalc = get_partiality(refl); if ( calc ) { @@ -426,8 +430,8 @@ static void partiality_histogram(cairo_t *cr, const struct image *images, } -static void scale_factor_histogram(cairo_t *cr, const struct image *images, - int n, const char *title) +static void scale_factor_histogram(cairo_t *cr, Crystal **crystals, int n, + const char *title) { int f_max; int i, b; @@ -451,8 +455,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, osf_max = 0.0; for ( i=0; i<n; i++ ) { - double osf = images[i].osf; - if ( images[i].pr_dud ) continue; + double osf = crystal_get_osf(crystals[i]); + if ( crystal_get_user_flag(crystals[i]) ) continue; if ( osf > osf_max ) osf_max = osf; } osf_max = ceil(osf_max+osf_max/10000.0); @@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, for ( i=0; i<n; i++ ) { - double osf = images[i].osf; + double osf = crystal_get_osf(crystals[i]); - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(crystals[i]) ) continue; for ( b=0; b<nbins; b++ ) { if ( (osf >= osf_low[b]) @@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, } -static void intensity_histogram(cairo_t *cr, const struct image *images, - int n, RefList *full, +static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n, + RefList *full, signed int h, signed int k, signed int l) { int f_max; @@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images, for ( i=0; i<n; i++ ) { double osf; + Crystal *cryst = crystals[i]; - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(cryst) ) continue; - osf = images[i].osf; + osf = crystal_get_osf(cryst); - for ( f = find_refl(images[i].reflections, h, k, l); + for ( f = find_refl(crystal_get_reflections(cryst), h, k, l); f != NULL; f = next_found_refl(f) ) { @@ -616,12 +621,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images, for ( i=0; i<n; i++ ) { double osf; + Crystal *cryst = crystals[i]; - if ( images[i].pr_dud ) continue; + if ( crystal_get_user_flag(cryst) ) continue; - osf = images[i].osf; + osf = crystal_get_osf(cryst); - for ( f = find_refl(images[i].reflections, h, k, l); + for ( f = find_refl(crystal_get_reflections(cryst), h, k, l); f != NULL; f = next_found_refl(f) ) { @@ -771,7 +777,7 @@ static void find_most_sampled_reflections(RefList *list, int n, signed int *h, -SRContext *sr_titlepage(struct image *images, int n, +SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename, const char *stream_filename, const char *cmdline) { @@ -805,7 +811,7 @@ SRContext *sr_titlepage(struct image *images, int n, } -void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, +void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n, RefList *full) { int i; @@ -822,7 +828,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 480.0, 350.0); - scale_factor_histogram(sr->cr, images, n, + scale_factor_histogram(sr->cr, crystals, n, "Distribution of overall scale factors"); cairo_restore(sr->cr); @@ -830,7 +836,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 70.0, 330.0); - partiality_graph(sr->cr, images, n, full); + partiality_graph(sr->cr, crystals, n, full); cairo_save(sr->cr); cairo_move_to(sr->cr, 0.0, 0.0); @@ -841,7 +847,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_stroke(sr->cr); cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 0.0, -150.0); - partiality_histogram(sr->cr, images, n, full, 1, 0); + partiality_histogram(sr->cr, crystals, n, full, 1, 0); cairo_restore(sr->cr); cairo_save(sr->cr); @@ -854,7 +860,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 230.0, 200.0); cairo_rotate(sr->cr, -M_PI_2); - partiality_histogram(sr->cr, images, n, full, 0, 1); + partiality_histogram(sr->cr, crystals, n, full, 0, 1); cairo_restore(sr->cr); cairo_restore(sr->cr); @@ -873,7 +879,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y); - intensity_histogram(sr->cr, images, n, full, + intensity_histogram(sr->cr, crystals, n, full, sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]); cairo_restore(sr->cr); diff --git a/src/scaling-report.h b/src/scaling-report.h index b9ac3fb7..5b153377 100644 --- a/src/scaling-report.h +++ b/src/scaling-report.h @@ -40,25 +40,25 @@ typedef struct _srcontext SRContext; /* Opaque */ #ifdef HAVE_CAIRO -extern SRContext *sr_titlepage(struct image *images, int n, +extern SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename, const char *stream_filename, const char *cmdline); -extern void sr_iteration(SRContext *sr, int iteration, struct image *images, - int n, RefList *full); +extern void sr_iteration(SRContext *sr, int iteration, + Crystal **crystals, int n, RefList *full); extern void sr_finish(SRContext *sr); #else -SRContext *sr_titlepage(struct image *images, int n, const char *filename, +SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename, const char *stream_filename, const char *cmdline) { return NULL; } -void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, +void sr_iteration(SRContext *sr, int iteration, Crystal **crystals, int n, RefList *full) { } diff --git a/tests/first_merge_check b/tests/first_merge_check index 686222cb..d649e284 100755 --- a/tests/first_merge_check +++ b/tests/first_merge_check @@ -1,34 +1,36 @@ #!/bin/sh cat > first_merge_check.stream << EOF -CrystFEL stream format 2.0 +CrystFEL stream format 2.1 Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 100.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 200.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- EOF diff --git a/tests/fourth_merge_check b/tests/fourth_merge_check index 619afc02..e2d2487f 100755 --- a/tests/fourth_merge_check +++ b/tests/fourth_merge_check @@ -1,20 +1,21 @@ #!/bin/sh cat > fourth_merge_check.stream << EOF -CrystFEL stream format 2.0 +CrystFEL stream format 2.1 Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 100.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- EOF diff --git a/tests/integration_check.c b/tests/integration_check.c index 7962e414..659eb7cc 100644 --- a/tests/integration_check.c +++ b/tests/integration_check.c @@ -205,7 +205,8 @@ int main(int argc, char *argv[]) image.height = 128; memset(image.data, 0, 128*128*sizeof(float)); - image.reflections=reflist_new(); + image.n_crystals = 0; + image.crystals = NULL; /* First check: no intensity -> no peak, or very low intensity */ r = integrate_peak(&image, 64, 64, &fsp, &ssp, &intensity, &sigma, diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index b39dfd05..9d713e0b 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -122,41 +122,98 @@ static void shift_parameter(struct image *image, int k, double shift) } -static void calc_either_side(struct image *image, double incr_val, +static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) +{ + Crystal *cr_new; + double r; + UnitCell *cell; + + cr_new = crystal_new(); + if ( cr_new == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return NULL; + } + + crystal_set_image(cr_new, crystal_get_image(cr)); + r = crystal_get_profile_radius(cr_new); + + switch ( refine ) { + + case REF_ASX : + case REF_ASY : + case REF_ASZ : + case REF_BSX : + case REF_BSY : + case REF_BSZ : + case REF_CSX : + case REF_CSY : + case REF_CSZ : + cell = new_shifted_cell(crystal_get_cell(cr), refine, + -incr_val); + crystal_set_cell(cr_new, cell); + break; + + case REF_R : + cell = cell_new_from_cell(crystal_get_cell(cr)); + crystal_set_profile_radius(cr_new, r + incr_val); + break; + + default : + ERROR("Can't shift %i\n", refine); + break; + + } + + return cr_new; +} + +static void calc_either_side(Crystal *cr, double incr_val, int *valid, long double *vals[3], int refine) { RefList *compare; - UnitCell *cell; + struct image *image = crystal_get_image(cr); + + if ( (refine != REF_DIV) ) { + + /* Crystal properties */ - if ( (refine != REF_DIV) && (refine != REF_R) ) { + Crystal *cr_new; - cell = new_shifted_cell(image->indexed_cell, refine, -incr_val); - compare = find_intersections(image, cell); - scan_partialities(image->reflections, compare, valid, vals, 0); - cell_free(cell); + cr_new = new_shifted_crystal(cr, refine, -incr_val); + + compare = find_intersections(image, cr_new); + scan_partialities(crystal_get_reflections(cr), compare, valid, + vals, 0); + cell_free(crystal_get_cell(cr_new)); + crystal_free(cr_new); reflist_free(compare); - cell = new_shifted_cell(image->indexed_cell, refine, +incr_val); - compare = find_intersections(image, cell); - scan_partialities(image->reflections, compare, valid, vals, 2); - cell_free(cell); + cr_new = new_shifted_crystal(cr, refine, +incr_val); + + compare = find_intersections(image, cr_new); + scan_partialities(crystal_get_reflections(cr), compare, valid, + vals, 0); + cell_free(crystal_get_cell(cr_new)); + crystal_free(cr_new); reflist_free(compare); } else { + /* "Image" properties */ + struct image im_moved; im_moved = *image; shift_parameter(&im_moved, refine, -incr_val); - compare = find_intersections(&im_moved, im_moved.indexed_cell); - scan_partialities(im_moved.reflections, compare, + compare = find_intersections(&im_moved, cr); + scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 0); reflist_free(compare); im_moved = *image; shift_parameter(&im_moved, refine, +incr_val); - compare = find_intersections(&im_moved, im_moved.indexed_cell); - scan_partialities(im_moved.reflections, compare, + compare = find_intersections(&im_moved, cr); + scan_partialities(crystal_get_reflections(cr), compare, valid, vals, 2); reflist_free(compare); @@ -165,7 +222,7 @@ static void calc_either_side(struct image *image, double incr_val, -static int test_gradients(struct image *image, double incr_val, int refine, +static int test_gradients(Crystal *cr, double incr_val, int refine, const char *str) { Reflection *refl; @@ -175,11 +232,13 @@ static int test_gradients(struct image *image, double incr_val, int refine, int *valid; int nref; int n_acc, n_valid; + RefList *reflections; //FILE *fh; - image->reflections = find_intersections(image, image->indexed_cell); + reflections = find_intersections(crystal_get_image(cr), cr); + crystal_set_reflections(cr, reflections); - nref = num_reflections(image->reflections); + nref = num_reflections(reflections); if ( nref < 10 ) { ERROR("Too few reflections found. Failing test by default.\n"); return -1; @@ -200,16 +259,15 @@ static int test_gradients(struct image *image, double incr_val, int refine, } for ( i=0; i<nref; i++ ) valid[i] = 1; - scan_partialities(image->reflections, image->reflections, - valid, vals, 1); + scan_partialities(reflections, reflections, valid, vals, 1); - calc_either_side(image, incr_val, valid, vals, refine); + calc_either_side(cr, incr_val, valid, vals, refine); //fh = fopen("wrongness.dat", "a"); n_valid = nref; n_acc = 0; i = 0; - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(reflections, &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -231,8 +289,7 @@ static int test_gradients(struct image *image, double incr_val, int refine, grad2 = (vals[2][i] - vals[1][i]) / incr_val; grad = (grad1 + grad2) / 2.0; - cgrad = gradient(image, refine, refl, - image->profile_radius); + cgrad = gradient(cr, refine, refl); get_partial(refl, &r1, &r2, &p, &cl, &ch); @@ -289,6 +346,7 @@ int main(int argc, char *argv[]) double bx, by, bz; double cx, cy, cz; UnitCell *cell; + Crystal *cr; struct quaternion orientation; int i; int val; @@ -303,10 +361,17 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J(8000.0)); image.div = 1e-3; image.bw = 0.01; - image.m = 0.0; - image.profile_radius = 0.005e9; image.filename = malloc(256); + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 1; + } + crystal_set_mosaicity(cr, 0.0); + crystal_set_profile_radius(cr, 0.005e9); + crystal_set_image(cr, &image); + cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9, deg2rad(90.0), deg2rad(90.0), @@ -316,36 +381,39 @@ int main(int argc, char *argv[]) for ( i=0; i<1; i++ ) { + UnitCell *rot; + orientation = random_quaternion(); - image.indexed_cell = cell_rotate(cell, orientation); + rot = cell_rotate(cell, orientation); + crystal_set_cell(cr, rot); - cell_get_reciprocal(image.indexed_cell, + cell_get_reciprocal(rot, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); incr_val = incr_frac * image.div; - val += test_gradients(&image, incr_val, REF_DIV, "div"); + val += test_gradients(cr, incr_val, REF_DIV, "div"); incr_val = incr_frac * ax; - val += test_gradients(&image, incr_val, REF_ASX, "ax*"); + val += test_gradients(cr, incr_val, REF_ASX, "ax*"); incr_val = incr_frac * ay; - val += test_gradients(&image, incr_val, REF_ASY, "ay*"); + val += test_gradients(cr, incr_val, REF_ASY, "ay*"); incr_val = incr_frac * az; - val += test_gradients(&image, incr_val, REF_ASZ, "az*"); + val += test_gradients(cr, incr_val, REF_ASZ, "az*"); incr_val = incr_frac * bx; - val += test_gradients(&image, incr_val, REF_BSX, "bx*"); + val += test_gradients(cr, incr_val, REF_BSX, "bx*"); incr_val = incr_frac * by; - val += test_gradients(&image, incr_val, REF_BSY, "by*"); + val += test_gradients(cr, incr_val, REF_BSY, "by*"); incr_val = incr_frac * bz; - val += test_gradients(&image, incr_val, REF_BSZ, "bz*"); + val += test_gradients(cr, incr_val, REF_BSZ, "bz*"); incr_val = incr_frac * cx; - val += test_gradients(&image, incr_val, REF_CSX, "cx*"); + val += test_gradients(cr, incr_val, REF_CSX, "cx*"); incr_val = incr_frac * cy; - val += test_gradients(&image, incr_val, REF_CSY, "cy*"); + val += test_gradients(cr, incr_val, REF_CSY, "cy*"); incr_val = incr_frac * cz; - val += test_gradients(&image, incr_val, REF_CSZ, "cz*"); + val += test_gradients(cr, incr_val, REF_CSZ, "cz*"); } diff --git a/tests/second_merge_check b/tests/second_merge_check index af952030..aec5b120 100755 --- a/tests/second_merge_check +++ b/tests/second_merge_check @@ -1,34 +1,36 @@ #!/bin/sh cat > second_merge_check.stream << EOF -CrystFEL stream format 2.0 +CrystFEL stream format 2.1 Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 100.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px -1 0 0 200.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- EOF diff --git a/tests/third_merge_check b/tests/third_merge_check index f3cfb637..569d5ddf 100755 --- a/tests/third_merge_check +++ b/tests/third_merge_check @@ -1,48 +1,51 @@ #!/bin/sh cat > third_merge_check.stream << EOF -CrystFEL stream format 2.0 +CrystFEL stream format 2.1 Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 100.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 200.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- ----- Begin chunk ----- Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 -I0 = 1.0 (arbitrary units) -photon_energy_eV = 2000.0 Reflections measured after indexing h k l I phase sigma(I) counts fs/px ss/px 1 0 0 100.00 - 0.00 1 938.0 629.0 End of reflections +--- End crystal ----- End chunk ----- EOF |