116 lines
3.9 KiB
C++
116 lines
3.9 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#ifndef JFJOCH_EIGENREFINE_H
|
|
#define JFJOCH_EIGENREFINE_H
|
|
|
|
#include <vector>
|
|
#include <optional>
|
|
|
|
#include "../common/CrystalLattice.h"
|
|
#include <ffbidx/refine.h>
|
|
|
|
struct RefineParameters {
|
|
int64_t viable_cell_min_spots;
|
|
float dist_tolerance_vs_reference;
|
|
std::optional<UnitCell> reference_unit_cell;
|
|
float min_length_A;
|
|
float max_length_A;
|
|
float indexing_tolerance;
|
|
};
|
|
|
|
inline std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
|
|
size_t nspots,
|
|
Eigen::MatrixX3<float> &oCell,
|
|
Eigen::VectorX<float> &scores,
|
|
RefineParameters &p) {
|
|
using M3x = Eigen::MatrixX3<float>;
|
|
|
|
std::vector<CrystalLattice> ret;
|
|
|
|
Eigen::MatrixX3<float> spots(in_spots.size(), 3u);
|
|
|
|
for (int i = 0; i < in_spots.size(); i++) {
|
|
spots(i, 0u) = in_spots[i].x;
|
|
spots(i, 1u) = in_spots[i].y;
|
|
spots(i, 2u) = in_spots[i].z;
|
|
}
|
|
|
|
fast_feedback::refine::config_ifssr<float> cifssr{
|
|
.min_spots = static_cast<uint32_t>(p.viable_cell_min_spots)
|
|
};
|
|
|
|
fast_feedback::refine::indexer_ifssr<float>::refine(spots.topRows(nspots), oCell, scores, cifssr);
|
|
|
|
// Select cell that explains most spots
|
|
int64_t max_indexed_spot_count = 0;
|
|
float min_score = -1;
|
|
int64_t id = -1;
|
|
|
|
for (int i = 0; i < scores.size(); i++) {
|
|
Eigen::Vector3f row_norms = oCell.block(3u * i, 0u, 3u, 3u).rowwise().norm();
|
|
|
|
// Check for distance vs. reference unit cell
|
|
if (p.reference_unit_cell) {
|
|
// Compare edge lengths up to 5% deviation, permutation-invariant
|
|
std::array<float, 3> obs = {row_norms(0), row_norms(1), row_norms(2)};
|
|
std::array<float, 3> ref = {
|
|
static_cast<float>(p.reference_unit_cell->a),
|
|
static_cast<float>(p.reference_unit_cell->b),
|
|
static_cast<float>(p.reference_unit_cell->c)
|
|
};
|
|
std::sort(obs.begin(), obs.end());
|
|
std::sort(ref.begin(), ref.end());
|
|
|
|
bool lengths_ok = true;
|
|
for (int k = 0; k < 3; ++k) {
|
|
// Guard against zero/near-zero reference values
|
|
const float denom = std::max(ref[k], 1e-6f);
|
|
const float rel_dev = std::abs(obs[k] - ref[k]) / denom;
|
|
if (rel_dev > p.dist_tolerance_vs_reference) {
|
|
lengths_ok = false;
|
|
break;
|
|
}
|
|
}
|
|
if (!lengths_ok) continue;
|
|
} else {
|
|
// Don't include lattices with weird lengths (only if no reference unit cell)
|
|
if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A)
|
|
continue;
|
|
}
|
|
|
|
// Count spots that match indexing tolerance
|
|
auto cell = oCell.block(3u * i, 0u, 3u, 3u);
|
|
M3x resid = spots.topRows(nspots) * cell.transpose();
|
|
const M3x miller = round(resid.array());
|
|
resid -= miller;
|
|
|
|
int64_t indexed_spot_count = (resid.rowwise().norm().array() < p.indexing_tolerance).count();
|
|
if (indexed_spot_count > max_indexed_spot_count) {
|
|
max_indexed_spot_count = indexed_spot_count;
|
|
min_score = scores(i);
|
|
id = i;
|
|
} if (indexed_spot_count == max_indexed_spot_count && scores(i) < min_score) {
|
|
min_score = scores(i);
|
|
id = i;
|
|
}
|
|
}
|
|
|
|
if (id == -1)
|
|
return {};
|
|
|
|
auto cell = oCell.block(3u * id, 0u, 3u, 3u);
|
|
|
|
fast_feedback::refine::make_right_handed(cell);
|
|
|
|
return { CrystalLattice(
|
|
Coord(cell(0,0), cell(0,1), cell(0,2)),
|
|
Coord(cell(1,0), cell(1,1), cell(1,2)),
|
|
Coord(cell(2,0), cell(2,1), cell(2,2))
|
|
)};
|
|
}
|
|
|
|
|
|
|
|
#endif //JFJOCH_EIGENREFINE_H
|