diff --git a/common/CrystalLattice.cpp b/common/CrystalLattice.cpp index ac3583c50..7a79630ae 100644 --- a/common/CrystalLattice.cpp +++ b/common/CrystalLattice.cpp @@ -72,6 +72,13 @@ float CrystalLattice::CalcVolume() const { return vec[0] * cross_product; } +float CrystalLattice::VolumeFraction() const { + const float denom = vec[0].Length() * vec[1].Length() * vec[2].Length(); + if (!(denom > 0.0f)) + return 0.0f; + return std::fabs(CalcVolume()) / denom; +} + void CrystalLattice::Sort() { if (vec[0].Length() > vec[1].Length()) std::swap(vec[0], vec[1]); @@ -94,15 +101,27 @@ void CrystalLattice::FixHandedness() { } +// The reciprocal basis divides by the cell volume, so three coplanar rows make all three vectors +// infinite - and an infinite a* does not fail, it quietly predicts nothing and poisons every +// residual built from it. Say so instead. +void CrystalLattice::CheckHasReciprocal() const { + if (VolumeFraction() < MIN_BASIS_VOLUME_FRACTION) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Crystal lattice is coplanar and has no reciprocal cell"); +} + Coord CrystalLattice::Astar() const { + CheckHasReciprocal(); return (vec[1] % vec[2]) * (1.0f / CalcVolume()); } Coord CrystalLattice::Bstar() const { + CheckHasReciprocal(); return (vec[2] % vec[0]) * (1.0f / CalcVolume()); } Coord CrystalLattice::Cstar() const { + CheckHasReciprocal(); return (vec[0] % vec[1]) * (1.0f / CalcVolume()); } diff --git a/common/CrystalLattice.h b/common/CrystalLattice.h index cb18b179a..96a38e0f9 100644 --- a/common/CrystalLattice.h +++ b/common/CrystalLattice.h @@ -11,9 +11,17 @@ #include "Coord.h" #include "UnitCell.h" +// Below this a basis counts as coplanar. 0.02 is 1.1 deg off flat: an order of magnitude below the +// flattest basis a real indexing candidate has been seen to reach, and three orders above the point +// where float cell angles stop carrying even the SIGN of the metric determinant. +constexpr float MIN_BASIS_VOLUME_FRACTION = 0.02f; + class CrystalLattice { Coord vec[3]; void FixHandedness(); + // Throws if the basis is too flat to have a reciprocal cell. Everything that builds a lattice + // rejects one this flat, so reaching it means a degenerate basis got through somewhere. + void CheckHasReciprocal() const; public: void ReorderABEqual(); void ReorderMonoclinic(); @@ -23,6 +31,11 @@ public: CrystalLattice(float a, float b, float c, float alpha, float beta, float gamma); CrystalLattice(const Coord &a, const Coord &b, const Coord &c); [[nodiscard]] float CalcVolume() const; + // |V| / (|a| |b| |c|) - the cell volume made scale-free: 1 for an orthogonal basis, 0 for three + // coplanar rows. This is the test for a degenerate basis, not CalcVolume(): three nearly-coplanar + // 300 A rows still enclose hundreds of A^3, so an absolute volume cannot see them, and a basis + // that flat has no usable reciprocal cell (1/V blows up) whatever its size. + [[nodiscard]] float VolumeFraction() const; [[nodiscard]] const Coord &Vec0() const; [[nodiscard]] const Coord &Vec1() const; [[nodiscard]] const Coord &Vec2() const; diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index ce5ad34d6..0976338e1 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -60,7 +60,7 @@ IndexingSettings &IndexingSettings::FFT_MaxUnitCell_A(float input) { IndexingSettings &IndexingSettings::FFT_MinUnitCell_A(float input) { check_finite("FFT indexing min unit cell (A)", input); - check_min("FFT indexing min unit cell (A)", input, 5); + check_min("FFT indexing min unit cell (A)", input, fft_min_unit_cell_limit_A); check_max("FFT indexing min unit cell (A)", input, 40); fft_min_unit_cell_A = input; return *this; diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 34e2d864b..eb139d568 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -33,6 +33,10 @@ public: // this rather than let the setter throw - a rescue that recovers an implausible axis must not // take the whole run down with it. static constexpr float fft_max_unit_cell_limit_A = 1200.0; + // The other end of the same search: any candidate whose reduced cell has an axis shorter than + // fft_min_unit_cell_A is discarded, and this is the lowest floor the setter accepts. A caller + // lowering the floor to reach a given cell clamps to it rather than let the setter throw. + static constexpr float fft_min_unit_cell_limit_A = 5.0; private: int64_t indexing_threads = 4; // Threads splitting the candidate-cell refinement WITHIN one indexer call. 1 (the default) is the diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index effba91c6..b6e72e6ad 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,9 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* A candidate cell whose three rows are coplanar is rejected before refinement, instead of producing a not-a-number Jacobian and several hundred lines of solver output. +* When the directions an FFT search shortlists all lie in one plane, one further transform is spent along the plane normal, which is where the missing row must be - so a crystal with a very long axis can index. +* `--fft-min-unit-cell` sets the shortest cell axis the FFT will consider, and giving a cell with `-C` lowers it automatically, so a small-molecule cell is no longer forced onto a supercell. * The lattice search rejects a Bravais character whose implied cell angle is geometrically impossible, instead of matching it because the comparison against a not-a-number is false. * The lattice search retries its Niggli type-boundary supplement for each of the three cell angles, not only for beta. * Systematic-absence evidence for a screw axis is scored per axial row rather than pooled, so an unmeasured row no longer vetoes a confirmed one. diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 7ce6079e8..53e7d15eb 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -137,6 +137,16 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, const float tolerance, const int num_threads) { try { + // A coplanar basis has no reciprocal cell: 1/V is infinite, every predicted reciprocal vector + // comes out NaN, and Ceres fails on the very first evaluation - after dumping the offending + // block to stderr. There is nothing for the refinement to recover here, so refuse the lattice + // before the problem is built rather than let the solver discover it. The check has to be on + // the vectors: this close to flat, float cell angles no longer carry even the SIGN of the + // metric determinant, and the triclinic branch of XtalResidual then clamps c into the a-b + // plane and divides by the zero volume that makes. + if (data.latt.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION) + return false; + Coord vec0 = data.latt.Vec0(); Coord vec1 = data.latt.Vec1(); Coord vec2 = data.latt.Vec2(); @@ -473,6 +483,11 @@ bool XtalOptimizerRotationOnly(XtalOptimizerData &data, const std::vector &spots, const float tolerance) { try { + // Same refusal as XtalOptimizerInternal: the residual here is built from Astar/Bstar/Cstar, + // which divide by the cell volume, so a coplanar basis makes every one of them infinite. + if (data.latt.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION) + return false; + // Parameter: angle-axis for the extra rotation. Identity == {0,0,0}. double rot_aa[3] = {0.0, 0.0, 0.0}; diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index 55098827b..5e49aa1df 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -12,6 +12,7 @@ #include "gemmi/symmetry.hpp" #include "../../common/JFJochException.h" +#include "../../common/CrystalLattice.h" // Rodrigues rotation with everything that depends only on the angle-axis worked out once. This is // ceres::AngleAxisRotatePoint term for term - so the rotated point is bit-identical to it - but that @@ -281,7 +282,14 @@ struct XtalResidual { cxa = c_unrot.cross(a_unrot); axb = a_unrot.cross(b_unrot); - invV = C(1) / a_unrot.dot(bxc); + // The clamp on cz above puts c in the a-b plane when the refined angles stop describing a + // real cell, and the triple product is then zero: 1/V is infinite and every residual built + // from it is NaN, which Ceres reports as a failed evaluation after dumping the block to + // stderr. Hold the volume at the same fraction of |a||b||c| that the candidate filter calls + // coplanar, so the step stays finite and the solver can walk back out of it. + const C V = a_unrot.dot(bxc); + const C V_min = C(MIN_BASIS_VOLUME_FRACTION) * L0 * L1 * L2; + invV = C(1) / ((V > V_min) ? V : V_min); } // h a* + k b* + l c* in the unrotated frame, from a basis ReciprocalBasis has already built. diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 8596507af..e994c6a94 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -124,6 +124,13 @@ std::vector FFTIndexer::ReduceResults(const std::vector & gamma < min_angle_deg || gamma > max_angle_deg) continue; + // Three coplanar rows are invisible to the two tests above: they give perfectly + // ordinary lengths and angles - 30 to 150 deg admits any flat combination - and + // only the volume betrays them. Such a cell has no reciprocal basis at all + // (1/V -> infinity), so it cannot be refined, only crashed into. + if (reduced.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION) + continue; + candidates.emplace_back(std::move(reduced)); } } @@ -151,6 +158,8 @@ std::vector FFTIndexer::ReduceResults(const std::vector & uc.beta < min_angle_deg || uc.beta > max_angle_deg || uc.gamma < min_angle_deg || uc.gamma > max_angle_deg) continue; + if (reduced.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION) + continue; bool duplicate = false; for (const auto &c : candidate_cells) @@ -340,6 +349,97 @@ std::vector FFTIndexer::ReduceAndRefine(const std::vector return Refine(coord, nspots, oCell, scores, parameters); } +// Vector perpendicular to a shortlist that lies in a single plane. Three vectors of one plane can +// never be reduced to a cell - every triple ReduceResults forms from them is degenerate - so a +// shortlist like that means the row that closes the cell was not among the search directions at all. +// The scatter matrix's smallest eigenvector is that plane's normal, and its eigenvalue ratio says how +// flat the set is: lambda_min/lambda_max = 2 . The threshold below is a +// 3.5 deg rms out-of-plane spread, which sits in a wide empty gap - measured over the accumulated +// clouds of ten rotation first passes, a degenerate shortlist scores 2e-5 to 3e-4 (0.2-0.8 deg) and +// every non-degenerate one 0.026 or more (7.7 deg and up). +std::optional FFTIndexer::DegeneratePlaneNormal(const std::vector &filtered) const { + // Below five rows the test says nothing: three directions drawn at random pass it 23% of the time + // and four 4.2%, against 0.7% at five and under 0.01% at eight (measured, 20000 draws each). A net + // that really is planar contributes far more rows than five - the crystal here gives sixteen. + if (filtered.size() < 5) + return {}; + + Eigen::Matrix3f scatter = Eigen::Matrix3f::Zero(); + for (const auto &v : filtered) { + const Coord n = v.Normalize(); + const Eigen::Vector3f e(n.x, n.y, n.z); + scatter += e * e.transpose(); + } + + Eigen::SelfAdjointEigenSolver es(scatter); + if (es.info() != Eigen::Success) + return {}; + + constexpr float MAX_OUT_OF_PLANE = 0.0076f; // 2* at eps_rms = 3.5 deg + if (es.eigenvalues()(0) > MAX_OUT_OF_PLANE * es.eigenvalues()(2)) + return {}; + + const auto v = es.eigenvectors().col(0); + return Coord(v(0), v(1), v(2)); +} + +// One more transform, with every search direction inside a cap about `axis`. The peak of a length-L +// row stays coherent only within roughly 10.7/(L * sigma_perp) degrees of its own direction (sigma_perp +// = the spread of the accumulated cloud across that direction), which for a long axis is finer than the +// ~0.43 deg the half-sphere grid leaves between a direction and its nearest node. Spending the same +// number of directions on a 3 deg cap samples it to ~0.016 deg instead, so the row is found where the +// global pass could only miss it. The count is unchanged, so the transform and its plan are untouched +// and this costs exactly one extra pass. +std::optional FFTIndexer::SearchCap(const std::vector &coord, size_t nspots, + const Coord &axis, float half_angle_deg) { + const Coord u = axis.Normalize(); + Coord seed = std::fabs(u.x) < 0.9f ? Coord(1.0f, 0.0f, 0.0f) : Coord(0.0f, 1.0f, 0.0f); + const Coord e1 = (u % seed).Normalize(); + const Coord e2 = u % e1; + + const double cos_half = std::cos(half_angle_deg * PI / 180.0); + const double golden_angle = 2.0 * PI / ((1.0 + std::sqrt(5.0)) / 2.0); + + std::vector saved; + saved.swap(direction_vectors); + direction_vectors.reserve(saved.size()); + for (size_t i = 0; i < saved.size(); i++) { + // Same spiral as SetupDirectionVectors, with z restricted to the cap instead of the half-sphere. + const double z = 1.0 - (1.0 - cos_half) * static_cast(i) / static_cast(saved.size() - 1); + const double r = std::sqrt(std::max(0.0, 1.0 - z * z)); + const double theta = golden_angle * static_cast(i); + direction_vectors.emplace_back(u * static_cast(z) + + e1 * static_cast(r * std::cos(theta)) + + e2 * static_cast(r * std::sin(theta))); + } + DirectionsChanged(); + + // The search borrows a member the object keeps using afterwards, and an indexer outlives a failed + // transform - the broker's pool logs the error and hands the same indexer the next image - so the + // grid has to go back even when ExecuteFFT throws. + try { + ExecuteFFT(coord, nspots); + } catch (...) { + direction_vectors.swap(saved); + DirectionsChanged(); + throw; + } + + int best = 0; + for (int i = 1; i < nDirections; i++) + if (result_fft[i].magnitude > result_fft[best].magnitude) + best = i; + const FFTResult peak = result_fft[best]; + const Coord found = direction_vectors.at(peak.direction) * peak.length; + + direction_vectors.swap(saved); + DirectionsChanged(); + + if (peak.length < min_length_A) + return {}; + return found; +} + std::vector FFTIndexer::RunInternal(const std::vector &coord, size_t nspots) { if (nspots > coord.size()) nspots = coord.size(); @@ -363,7 +463,15 @@ std::vector FFTIndexer::RunInternal(const std::vector &co // standard candidates (the widened pass never runs). const float frac = lattices.empty() ? 0.0f : IndexedFraction(lattices.front(), coord, nspots); if (frac < 0.5f) { - for (auto &w : ReduceAndRefine(coord, nspots, FilterFFTResults(60), true)) { + auto filtered = FilterFFTResults(60); + // A shortlist confined to one plane cannot close a cell whatever is done with it, and the row + // that would close it is the one perpendicular to that plane. Look for it there, at an angular + // resolution the global grid does not have. + if (const auto normal = DegeneratePlaneNormal(filtered)) + if (const auto v = SearchCap(coord, nspots, *normal, 3.0f)) + filtered.push_back(*v); + + for (auto &w : ReduceAndRefine(coord, nspots, filtered, true)) { bool duplicate = false; for (const auto &l : lattices) if (l.GetUnitCell().is_close(w.GetUnitCell(), 0.02f, 1.0f)) { duplicate = true; break; } diff --git a/image_analysis/indexing/FFTIndexer.h b/image_analysis/indexing/FFTIndexer.h index b5061e3d4..dd0ab531b 100644 --- a/image_analysis/indexing/FFTIndexer.h +++ b/image_analysis/indexing/FFTIndexer.h @@ -42,8 +42,17 @@ protected: const std::vector &filtered, bool widen); float IndexedFraction(const CrystalLattice &latt, const std::vector &coord, size_t nspots) const; void SetupDirectionVectors(); + // Replace the search directions with a small cap around `axis` (same count, so the transform and + // its plan are unchanged) and return the strongest peak found in it. + std::optional SearchCap(const std::vector &coord, size_t nspots, + const Coord &axis, float half_angle_deg); + // Direction perpendicular to every vector of `filtered`, when they all lie in one plane and so + // cannot close a cell; empty otherwise. + std::optional DegeneratePlaneNormal(const std::vector &filtered) const; virtual void ExecuteFFT(const std::vector &coord, size_t nspots) = 0; + // Called after direction_vectors is rewritten, for implementations that keep a copy of it. + virtual void DirectionsChanged() {} public: explicit FFTIndexer(const IndexingSettings& settings); diff --git a/image_analysis/indexing/FFTIndexerGPU.cu b/image_analysis/indexing/FFTIndexerGPU.cu index 22b0b1400..343446845 100644 --- a/image_analysis/indexing/FFTIndexerGPU.cu +++ b/image_analysis/indexing/FFTIndexerGPU.cu @@ -175,19 +175,7 @@ FFTIndexerGPU::FFTIndexerGPU(const IndexingSettings &settings) d_dir_y = CudaDevicePtr(nDirections); d_dir_z = CudaDevicePtr(nDirections); - std::vector dir_x(nDirections); - std::vector dir_y(nDirections); - std::vector dir_z(nDirections); - - for (int i = 0; i < nDirections; i++) { - dir_x[i] = direction_vectors.at(i).x; - dir_y[i] = direction_vectors.at(i).y; - dir_z[i] = direction_vectors.at(i).z; - } - - cudaMemcpy(d_dir_x, dir_x.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); - cudaMemcpy(d_dir_y, dir_y.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); - cudaMemcpy(d_dir_z, dir_z.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); + DirectionsChanged(); int n[1] = {static_cast(histogram_size)}; // Size of the FFT along a single dimension @@ -197,6 +185,18 @@ FFTIndexerGPU::FFTIndexerGPU(const IndexingSettings &settings) } +void FFTIndexerGPU::DirectionsChanged() { + std::vector dir_x(nDirections), dir_y(nDirections), dir_z(nDirections); + for (int i = 0; i < nDirections; i++) { + dir_x[i] = direction_vectors.at(i).x; + dir_y[i] = direction_vectors.at(i).y; + dir_z[i] = direction_vectors.at(i).z; + } + cudaMemcpy(d_dir_x, dir_x.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(d_dir_y, dir_y.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); + cudaMemcpy(d_dir_z, dir_z.data(), nDirections * sizeof(float), cudaMemcpyHostToDevice); +} + void FFTIndexerGPU::ExecuteFFT(const std::vector &coord, size_t nspots) { int l_blockDim = 128; int l_gridDim = (direction_vectors.size() + l_blockDim - 1) / l_blockDim; diff --git a/image_analysis/indexing/FFTIndexerGPU.h b/image_analysis/indexing/FFTIndexerGPU.h index 0a1bdbbac..4f4938545 100644 --- a/image_analysis/indexing/FFTIndexerGPU.h +++ b/image_analysis/indexing/FFTIndexerGPU.h @@ -40,6 +40,7 @@ class FFTIndexerGPU : public FFTIndexer { CudaStream stream; void ExecuteFFT(const std::vector &coord, size_t nspots) override; + void DirectionsChanged() override; public: explicit FFTIndexerGPU(const IndexingSettings& settings); FFTIndexerGPU(const FFTIndexerGPU &i) = delete; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 5e0b63f3d..c10e62663 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -375,6 +375,9 @@ static void LogXDSGeometry(const DiffractionExperiment &experiment) { // with the margin DIALS uses when a cell is known (max_cell = 1.3 * longest axis, see // dials/algorithms/indexing/indexer.py find_max_cell). // +// The floor at the short end is adjusted here too, for the same reason and against the same cell - +// see the second half of the body. +// // ROTATION ONLY. Stills fire the FFT once per image across every worker and routinely run with a // known cell, so widening the transform there would cost the whole serial run for nothing; the // caller applies this inside its rotation branch. @@ -397,6 +400,21 @@ static void RaiseFFTBoundForKnownCell(IndexingSettings &settings, const Diffract Logger("Rugnux").Info("FFT search bound raised to {:.0f} A to cover the given cell " "(longest axis {:.1f} A)", needed, longest); } + + // The same argument at the other end of the search. Every candidate whose reduced cell has an + // axis below fft_min_unit_cell_A is discarded, and the 10 A default is shorter than any protein + // axis but longer than a small molecule's: a given cell with an 8 A axis was thrown away before + // it could be compared against the reference it was supposed to match. Lower the floor to reach + // the cell that was asked for, with the reciprocal of the margin used above. + constexpr float KNOWN_CELL_FLOOR_MARGIN = 1.0f / KNOWN_CELL_MARGIN; + const float shortest = std::min({cell->a, cell->b, cell->c}); + const float floor_A = std::max(KNOWN_CELL_FLOOR_MARGIN * shortest, + IndexingSettings::fft_min_unit_cell_limit_A); + if (floor_A < settings.GetFFT_MinUnitCell_A()) { + settings.FFT_MinUnitCell_A(floor_A); + Logger("Rugnux").Info("FFT shortest-axis floor lowered to {:.1f} A to cover the given cell " + "(shortest axis {:.1f} A)", floor_A, shortest); + } } void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) { diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 5918604fa..20a916b99 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -133,6 +133,7 @@ void print_usage() { std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; std::cout << " -S, --space-group Space group number (92) or symbol (P43212) - for indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; + std::cout << " --fft-min-unit-cell Shortest cell axis the FFT search accepts, in A (default: 10). A candidate with a shorter axis is discarded, so a crystal below the default cannot be indexed at all; -C lowers it on its own to cover the cell given" << std::endl; std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|flex); flex tries all three per image and keeps whichever indexes the most spots (alias: multi)" << std::endl; std::cout << " --refine-geometry[=N|off] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200) then re-indexes (lifts weak-stills indexing). Default ON for stills when a reference cell is given (-C / reference MTZ); =off disables" << std::endl; std::cout << " --index-ice-rings[=on|off] Index on the spots flagged as sitting on an ice ring too, instead of setting them aside (default: off; no effect without --detect-ice-rings)" << std::endl; @@ -285,6 +286,7 @@ enum { OPT_ROT1, OPT_ROT2, OPT_ROT3, + OPT_FFT_MIN_UNIT_CELL, OPT_POLARIZATION }; @@ -336,6 +338,7 @@ static option long_options[] = { {"rot1", required_argument, nullptr, OPT_ROT1}, {"rot2", required_argument, nullptr, OPT_ROT2}, {"rot3", required_argument, nullptr, OPT_ROT3}, + {"fft-min-unit-cell", required_argument, nullptr, OPT_FFT_MIN_UNIT_CELL}, {"polarization", required_argument, nullptr, OPT_POLARIZATION}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, {"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE}, @@ -660,6 +663,7 @@ static int RunRugnux(int argc, char **argv) { bool index_ice_rings = false; // --index-ice-rings[=on|off]; index on the ice-band spots too std::optional ice_min_score_arg; // --ice-min-score: ice-presence gate on the measured score std::optional ice_min_spot_ratio_arg; // --ice-min-spot-ratio: the same gate on the spot channel + std::optional fft_min_unit_cell_A; // --fft-min-unit-cell: shortest axis the FFT search accepts std::optional min_q, max_q, q_spacing; // azimuthal integration range / -q spacing (1/A) std::optional azimuthal_bins; // --azimuthal-bins std::optional polarization_correction; // --polarization-correction (azimuthal integration) @@ -1235,6 +1239,8 @@ static int RunRugnux(int argc, char **argv) { case OPT_ROT1: rot1_rad = parse_float_arg(optarg, "--rot1", logger); break; case OPT_ROT2: rot2_rad = parse_float_arg(optarg, "--rot2", logger); break; case OPT_ROT3: rot3_rad = parse_float_arg(optarg, "--rot3", logger); break; + case OPT_FFT_MIN_UNIT_CELL: + fft_min_unit_cell_A = parse_float_arg(optarg, "--fft-min-unit-cell", logger); break; case OPT_POLARIZATION: polarization_factor = parse_float_arg(optarg, "--polarization", logger); break; case OPT_SCALING_ITERATIONS: scaling_iter = atoi(optarg); @@ -2083,6 +2089,7 @@ static int RunRugnux(int argc, char **argv) { indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value()); indexing_settings.GeomRefinementAlgorithm(refinement_algorithm); indexing_settings.IndexIceRings(index_ice_rings); + if (fft_min_unit_cell_A) indexing_settings.FFT_MinUnitCell_A(fft_min_unit_cell_A.value()); experiment.ImportIndexingSettings(indexing_settings); // --detect-ice-rings[=on|off] overrides the value carried in from the dataset (HDF5MetadataSource diff --git a/tests/CrystalLatticeTest.cpp b/tests/CrystalLatticeTest.cpp index d99d57577..decb29f22 100644 --- a/tests/CrystalLatticeTest.cpp +++ b/tests/CrystalLatticeTest.cpp @@ -113,6 +113,26 @@ TEST_CASE("CrystalLattice_Volume") { REQUIRE(l2.CalcVolume() == Catch::Approx(50 * 60 * 80 * sin120)); } +TEST_CASE("CrystalLattice_VolumeFraction") { + CrystalLattice ortho(50, 60, 80, 90, 90, 90); + REQUIRE(ortho.VolumeFraction() == Catch::Approx(1.0)); + + // All three angles at 60 deg is about as oblique as a reduced cell gets, and is still an order + // of magnitude clear of MIN_BASIS_VOLUME_FRACTION. + CrystalLattice oblique(50, 50, 50, 60, 60, 60); + REQUIRE(oblique.VolumeFraction() == Catch::Approx(std::sqrt(0.5)).margin(1e-4)); + REQUIRE(oblique.VolumeFraction() > MIN_BASIS_VOLUME_FRACTION); + + // Three coplanar rows: ordinary lengths, ordinary angles, no volume. This is what an indexing + // candidate built from directions that all lie in one dense reciprocal plane looks like, and + // what nothing downstream can refine - 1/V is not finite. + CrystalLattice flat(Coord(50, 0, 0), Coord(0, 60, 0), Coord(30, 40, 0)); + REQUIRE(flat.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION); + + CrystalLattice degenerate(Coord(0, 0, 0), Coord(0, 60, 0), Coord(0, 0, 80)); + REQUIRE(degenerate.VolumeFraction() == 0.0f); +} + TEST_CASE("CrystalLattice_Recip") { CrystalLattice l(50,60,80, 90, 90, 90); REQUIRE(l.Astar().Length() == Catch::Approx(1/50.0));