diff --git a/docs/RUGNUX_CLI.md b/docs/RUGNUX_CLI.md index a53e0c49b..1e25bf517 100644 --- a/docs/RUGNUX_CLI.md +++ b/docs/RUGNUX_CLI.md @@ -256,7 +256,7 @@ For very long axes the direction grid's angular resolution binds well below 500 | `--force-still` | Treat a rotation (goniometer) dataset as independent stills instead of rotation | | `-X, --indexing-algorithm ` | `Auto` (default) \| `FFT` \| `cuFFT` \| `FFTW` \| `FFBIDX` \| `None`, case-insensitive. `FFT` is the de-novo FFT search, on the GPU (cuFFT) where there is one and `--no-gpu` is not given, otherwise on the CPU (FFTW); `cuFFT` and `FFTW` name one implementation each; `cuFFT` needs a GPU and is refused without one. `FFBIDX` (known cell, `-C`) runs on the GPU where there is one and otherwise on the CPU. `Auto` takes `FFBIDX` on stills with a known cell, and `FFT` otherwise; rotation data always use the FFT search | | `-F` | Short for `-X FFT`, kept for old command lines | -| `-C, --unit-cell ` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it | +| `-C, --unit-cell ` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it. A candidate that is an integer supercell of the cell given (a twin's coincidence lattice, a row taken at a multiple of its length) is stepped down to it | | `--fft-min-unit-cell ` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but a de-novo rotation run also runs a **second first-pass hypothesis with the floor lowered to 5 Å by default**, adopted when the standing cell turns out to be an integer supercell of what it finds, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given (and the second hypothesis then stays out of the way). On stills, where neither applies, a crystal below the floor cannot be indexed unless this is lowered | | `--min-indexed-spots ` | Spots a frame must have on the lattice before it counts as indexed (default: 6, minimum 4). It sets the reported indexing rate and the count the rotation first pass ranks candidate lattices by; whether the run has a lattice at all is decided on the pooled spots instead ([§4.1](CPU_DATA_ANALYSIS_INDEXING.md#41-indexed-spot-decision-inlier-test)), and integration is not gated by it | | `-S, --space-group ` | Space group number (`96`) or Hermann-Mauguin symbol (`P43212`) — for indexing and scaling. On serial stills whose group admits two indexings, it also lets the crystals be put on one from the data ([the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)) | diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index ea906ae12..918cef72e 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -7,6 +7,7 @@ #include #include #include +#include namespace { struct config_ifssr final { @@ -142,6 +143,106 @@ namespace { cells.block(3u * j, 0u, 3u, 3u) = cell.transpose(); } } + + // The reference test, on a candidate's three rows: sorted lengths within the relative tolerance and + // sorted angles within REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG of any of the reference cells. Each + // angle is folded to its acute complement (min(x,180-x)) so the obtuse/acute setting choice is + // irrelevant. Guards against a right-edges/wrong-angle cell (a pseudo-symmetric near-metric, e.g. a + // monoclinic beta refined to the wrong value) passing on lengths. + bool MatchesReference(const Eigen::Matrix3f &cell_rows, const std::vector &reference_cells, + float dist_tolerance) { + std::array obs = {cell_rows.row(0).norm(), cell_rows.row(1).norm(), cell_rows.row(2).norm()}; + std::sort(obs.begin(), obs.end()); + + auto fold = [](float deg) { return std::min(deg, 180.0f - deg); }; + auto row_angle = [&](int i, int j) { + return std::acos(std::clamp(cell_rows.row(i).normalized().dot(cell_rows.row(j).normalized()), + -1.0f, 1.0f)) * 180.0f / PI; + }; + std::array obs_ang = {fold(row_angle(1, 2)), fold(row_angle(0, 2)), fold(row_angle(0, 1))}; + std::sort(obs_ang.begin(), obs_ang.end()); + + auto matches = [&](const UnitCell &reference) { + std::array ref = {reference.a, reference.b, reference.c}; + std::sort(ref.begin(), ref.end()); + for (int k = 0; k < 3; ++k) { + const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON); + if (std::abs(obs[k] - ref[k]) / denom > dist_tolerance) + return false; + } + + std::array ref_ang = {fold(reference.alpha), fold(reference.beta), fold(reference.gamma)}; + std::sort(ref_ang.begin(), ref_ang.end()); + for (int k = 0; k < 3; ++k) { + if (std::abs(obs_ang[k] - ref_ang[k]) > REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG) + return false; + } + return true; + }; + return std::any_of(reference_cells.begin(), reference_cells.end(), matches); + } + + // A candidate that is an integer SUPERCELL of the reference lattice: the FFT builds one from rows it + // reports at a multiple of their length, and on a twin by reticular merohedry the cell indexing both + // domains is their coincidence lattice, a supercell of either. Every node of the candidate is then a + // node of the reference lattice, so each reference basis vector is a combination of the candidate's + // rows with coefficients in (1/n)Z, n the volume ratio. Those combinations are enumerated, and of the + // triples that match the reference, span its volume and contain the candidate's lattice, the one + // that indexes the most spots is returned (on a twin, the stronger domain). + std::optional StepDownToReference(const Eigen::Matrix3f &cell_rows, + const std::vector &reference_cells, + float dist_tolerance, + const Eigen::Ref> &spots, + float indexing_tolerance) { + constexpr int MAX_INDEX = 10; + const float vol = std::abs(cell_rows.determinant()); + std::optional best; + int64_t best_count = -1; + for (const auto &ref: reference_cells) { + const float ref_vol = CrystalLattice(ref).CalcVolume(); + if (!(ref_vol > 0.0f)) + continue; + const float ratio = vol / ref_vol; + const int n = static_cast(std::lround(ratio)); + if (n < 2 || n > MAX_INDEX || std::abs(ratio - static_cast(n)) > 0.15f * static_cast(n)) + continue; + + const std::array lengths = {ref.a, ref.b, ref.c}; + std::array, 3> near; + for (int i = -n; i <= n; i++) + for (int j = -n; j <= n; j++) + for (int k = -n; k <= n; k++) { + const Eigen::RowVector3f v = (static_cast(i) * cell_rows.row(0) + + static_cast(j) * cell_rows.row(1) + + static_cast(k) * cell_rows.row(2)) + / static_cast(n); + for (int a = 0; a < 3; a++) + if (std::abs(v.norm() - lengths[a]) <= dist_tolerance * lengths[a]) + near[a].push_back(v); + } + + for (const auto &v0: near[0]) + for (const auto &v1: near[1]) + for (const auto &v2: near[2]) { + Eigen::Matrix3f sub; + sub << v0, v1, v2; + if (std::abs(std::abs(sub.determinant()) - ref_vol) > 0.15f * ref_vol) + continue; + const Eigen::Matrix3f m = cell_rows * sub.inverse(); + if ((m - m.array().round().matrix()).cwiseAbs().maxCoeff() > 0.1f) + continue; + if (!MatchesReference(sub, {ref}, dist_tolerance)) + continue; + int64_t count = 0; + ComputeIndexedMask(spots, sub.transpose(), indexing_tolerance, count); + if (count > best_count) { + best_count = count; + best = sub; + } + } + } + return best; + } } std::vector Refine(const std::vector &in_spots, @@ -216,41 +317,15 @@ std::vector Refine(const std::vector &in_spots, Eigen::Vector3f row_norms = cell_rows.rowwise().norm(); if (!reference_cells.empty()) { - std::array obs = {row_norms(0), row_norms(1), row_norms(2)}; - std::sort(obs.begin(), obs.end()); - - // The angles are compared as well as the lengths. Each is folded to its acute complement - // (min(x,180-x)) so the obtuse/acute setting choice is irrelevant, then the sorted - // triples are compared. Guards against a right-edges/wrong-angle cell (a pseudo-symmetric - // near-metric, e.g. a monoclinic beta refined to the wrong value) passing on lengths. - auto fold = [](float deg) { return std::min(deg, 180.0f - deg); }; - auto row_angle = [&](int i, int j) { - return std::acos(std::clamp(cell_rows.row(i).normalized().dot(cell_rows.row(j).normalized()), - -1.0f, 1.0f)) * 180.0f / PI; - }; - std::array obs_ang = {fold(row_angle(1, 2)), fold(row_angle(0, 2)), fold(row_angle(0, 1))}; - std::sort(obs_ang.begin(), obs_ang.end()); - - auto matches = [&](const UnitCell &reference) { - std::array ref = {reference.a, reference.b, reference.c}; - std::sort(ref.begin(), ref.end()); - for (int k = 0; k < 3; ++k) { - const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON); - if (std::abs(obs[k] - ref[k]) / denom > p.dist_tolerance_vs_reference) - return false; - } - - std::array ref_ang = {fold(reference.alpha), fold(reference.beta), fold(reference.gamma)}; - std::sort(ref_ang.begin(), ref_ang.end()); - for (int k = 0; k < 3; ++k) { - if (std::abs(obs_ang[k] - ref_ang[k]) > REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG) - return false; - } - return true; - }; - - if (std::none_of(reference_cells.begin(), reference_cells.end(), matches)) - continue; + if (!MatchesReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference)) { + const auto sub = StepDownToReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference, + spots.topRows(nspots), p.indexing_tolerance); + if (!sub) + continue; + cell_rows = *sub; + cell_cols = cell_rows.transpose(); + row_norms = cell_rows.rowwise().norm(); + } } else { if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A) continue; diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index 535d0dcbd..953b2e4ea 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -280,6 +280,46 @@ TEST_CASE("FFTIndexer_AxisLongerThanTheBound","[Indexing]") { CHECK(out[2] == Catch::Approx(360.0f).epsilon(0.01)); } +TEST_CASE("PostIndexingRefinement_SupercellStepsDownToReference","[Indexing]") { + // A candidate that is a 2x supercell of the given cell (one axis doubled, as the coincidence lattice + // of a twin or a row reported at its first multiple gives) comes back as the given cell, not refused. + const UnitCell uc(6.0, 8.0, 11.0, 90, 100, 90); + const CrystalLattice cl(uc); + + std::vector spots; + for (int h = -5; h <= 5; h++) + for (int k = -6; k <= 6; k++) + for (int l = -8; l <= 8; l++) + if (h != 0 || k != 0 || l != 0) + spots.push_back(h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar()); + + Eigen::MatrixX3 cells(3, 3); + const Coord a2 = cl.Vec0() * 2.0f, b = cl.Vec1(), c = cl.Vec2(); + cells << a2.x, a2.y, a2.z, b.x, b.y, b.z, c.x, c.y, c.z; + Eigen::VectorX scores(1); + scores(0) = 0.2f; + + RefineParameters p{ + .viable_cell_min_spots = 6, + .dist_tolerance_vs_reference = 0.05f, + .reference_unit_cell = uc, + .min_length_A = 2.5f, + .max_length_A = 500.0f, + .min_angle_deg = 30.0f, + .max_angle_deg = 150.0f, + .indexing_tolerance = 0.1f + }; + const auto out = Refine(spots, spots.size(), cells, scores, p); + REQUIRE(out.size() == 1); + CHECK(std::abs(out[0].CalcVolume()) == Catch::Approx(cl.CalcVolume()).epsilon(0.01)); + const auto uc_out = out[0].GetUnitCell(); + std::array lengths = {uc_out.a, uc_out.b, uc_out.c}; + std::sort(lengths.begin(), lengths.end()); + CHECK(lengths[0] == Catch::Approx(6.0f).epsilon(0.01)); + CHECK(lengths[1] == Catch::Approx(8.0f).epsilon(0.01)); + CHECK(lengths[2] == Catch::Approx(11.0f).epsilon(0.01)); +} + TEST_CASE("FFTIndexer_CentredReferenceCell","[Indexing]") { // A reference cell is quoted the way a deposit states it, which for a centred lattice is the // conventional cell - and the indexer's candidates are Niggli-reduced primitive cells, whose