mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 72 Wiki Activity

jfjoch_viewer: Calibrant can be freely chosen from options (LaB6, Ag behenate, current sample)

Browse Source
This commit is contained in:
leonarski_f
2025-10-29 14:49:14 +01:00
parent effe97a970
commit 549c0fade1
7 changed files with 109 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/geom_refinement/AssignSpotsToRings.cpp
+25 -14
View File
@@ -9,6 +9,8 @@
#include <algorithm>
#include <iostream>
#include "../../common/CrystalLattice.h"
FindCircleCenterResult FindCircleCenter(const std::vector<SpotToSave> &v, int64_t width, int64_t height, int64_t max_spots) {
if ((width <= 0) || (height <= 0) || (max_spots <= 0))
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid image size");
@@ -148,16 +150,21 @@ std::vector<RingClusters> AnalyzeClusters(const std::vector<float>& r, const std
return ret;
}
// Build unique sorted Qu = sqrt(h^2 + k^2 + l^2)
// This handles rings spacing for (any) crystal with cubic symmetry
std::vector<float> CalculateCubicXtalRings(float a, int hkl_max) {
std::vector<float> CalculateXtalRings(const UnitCell &cell, int hkl_max) {
CrystalLattice latt(cell);
Coord Astar = latt.Astar();
Coord Bstar = latt.Bstar();
Coord Cstar = latt.Cstar();
std::vector<float> u;
for (int h = 0; h <= hkl_max; h++) {
for (int k = 0; k <= hkl_max; k++) {
for (int l = 0; l <= hkl_max; l++) {
if (h == 0 && k == 0 && l == 0) continue;
float val = std::sqrt(float(h*h + k*k + l*l));
u.push_back(val);
auto p = Astar * h + Bstar * k + Cstar * l;
float Q = 2.0f * M_PI * p.Length();
u.push_back(Q);
}
}
}
@@ -165,15 +172,17 @@ std::vector<float> CalculateCubicXtalRings(float a, int hkl_max) {
// Deduplicate (since e.g. (100), (010), (001) all give sqrt(1))
u.erase(std::unique(u.begin(), u.end(),
[](float a, float b){ return std::fabs(a-b) < 1e-12; }),
[](float a, float b){ return std::fabs(a-b) < 1e-6; }),
u.end());
for (auto &iter: u)
iter = 2 * M_PI * iter / a;
return u;
}
std::vector<float> CalculateCubicXtalRings(float a, int hkl_max) {
return CalculateXtalRings(UnitCell(a,a,a,90,90,90), hkl_max);
}
float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A) {
float sin_theta = geom.GetWavelength_A() / (2 * d_A);
if (sin_theta < 0 || sin_theta > 1)
@@ -209,10 +218,12 @@ std::vector<RingClusters> GuessInitialGeometry(DiffractionGeometry &geom, const
return cluster_annot;
}
void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, float calibrant_a_A) {
auto cluster_annot = GuessInitialGeometry(geom, v, calibrant_a_A);
void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, const UnitCell &calibrant_uc) {
std::vector<float> ring_q = CalculateXtalRings(calibrant_uc);
if (ring_q.empty())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Unit cell problem");
auto cluster_annot = GuessInitialGeometry(geom, v, 2 * M_PI / ring_q[0]);
std::vector<float> ring_q = CalculateCubicXtalRings(calibrant_a_A);
std::vector<RingOptimizerInput> optimizer_input;
int ring_idx = 0;
@@ -246,8 +257,8 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v,
geom = optimizer.Run(optimizer_input);
}
void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, float calibrant_a_A) {
std::vector<float> q_ring = CalculateCubicXtalRings(calibrant_a_A);
void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, const UnitCell &calibrant) {
std::vector<float> q_ring = CalculateXtalRings(calibrant);
std::vector<RingOptimizerInput> optimizer_input;
image_analysis/geom_refinement/AssignSpotsToRings.h
+4 -3
View File
@@ -5,7 +5,7 @@
#define JFJOCH_ASSIGNSPOTSTORINGS_H
#include <vector>
#include <optional>
#include "../../common/UnitCell.h"
#include "../../common/SpotToSave.h"
#include "RingOptimizer.h"
@@ -29,11 +29,12 @@ FindCircleCenterResult FindCircleCenter(const std::vector<SpotToSave> &v,
std::vector<std::vector<int>> ClusterSpotsIntoRings(const std::vector<float>& r, float eps = 1.5, int minPts = 5);
std::vector<RingClusters> AnalyzeClusters(const std::vector<float>& r, const std::vector<std::vector<int>> &clusters);
std::vector<float> CalculateXtalRings(const UnitCell &cell, int hkl_max = 6);
std::vector<float> CalculateCubicXtalRings( float a, int hkl_max = 4);
float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A);
std::vector<RingClusters> GuessInitialGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, float calibrant_a_A);
void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, float calibrant_a_A);
void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, float calibrant_a_A);
void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, const UnitCell &calibrant);
void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, const UnitCell &calibrant);
#endif //JFJOCH_ASSIGNSPOTSTORINGS_H
tests/DetGeomCalibTest.cpp
+16 -4
View File
@@ -98,8 +98,8 @@ TEST_CASE("DetGeomCalib_AnalyzeClusters") {
REQUIRE(ret[1].R_obs == Catch::Approx(100.0f));
}
TEST_CASE("DetGeomCalib_build_u") {
auto ret = CalculateCubicXtalRings(2.0);
TEST_CASE("DetGeomCalib_CalculateXtalRings_cubic") {
auto ret = CalculateXtalRings(UnitCell(2.0, 2.0, 2.0, 90, 90, 90));
CHECK(ret[0] == Catch::Approx(2.0 * M_PI * 1.0 / 2.0));
CHECK(ret[1] == Catch::Approx(2.0 * M_PI * sqrt( 2.0 )/ 2.0));
CHECK(ret[2] == Catch::Approx(2.0 * M_PI * sqrt( 3.0 )/ 2.0));
@@ -108,6 +108,18 @@ TEST_CASE("DetGeomCalib_build_u") {
CHECK(ret[6] == Catch::Approx(2.0 * M_PI * sqrt( 8.0 )/ 2.0));
}
TEST_CASE("DetGeomCalib_CalculateXtalRings_one_long_axis") {
auto ret = CalculateXtalRings(UnitCell(50.0, 2.0, 2.0, 90, 90, 90));
CHECK(ret[0] == Catch::Approx(2.0 * M_PI * 1.0 / 50.0));
CHECK(ret[1] == Catch::Approx(2.0 * M_PI * 2.0 / 50.0));
CHECK(ret[2] == Catch::Approx(2.0 * M_PI * 3.0 / 50.0));
CHECK(ret[3] == Catch::Approx(2.0 * M_PI * 4.0 / 50.0));
CHECK(ret[4] == Catch::Approx(2.0 * M_PI * 5.0 / 50.0));
}
TEST_CASE("DetGeomCalib_GuessDetectorDistance") {
std::vector<SpotToSave> spots;
@@ -178,7 +190,7 @@ TEST_CASE("DetGeomCalib_GuessGeometry") {
DiffractionGeometry geom_out;
geom_out.Wavelength_A(1.0).DetectorDistance_mm(200.0);
GuessGeometry(geom_out, spots, LAB6_CELL_A);
GuessGeometry(geom_out, spots, UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90));
CHECK(fabsf(geom_out.GetBeamX_pxl() - geom.GetBeamX_pxl()) < 0.001f);
CHECK(fabsf(geom_out.GetBeamY_pxl() - geom.GetBeamY_pxl()) < 0.001f);
@@ -246,7 +258,7 @@ TEST_CASE("DetGeomCalib_OptimizeGeometry") {
geom_i.Wavelength_A(1.0).BeamX_pxl(995.0).BeamY_pxl(1277.0)
.DetectorDistance_mm(98).PoniRot1_rad(0.0975).PoniRot2_rad(0.055);
OptimizeGeometry(geom_i, spots, lab6_a);
OptimizeGeometry(geom_i, spots, UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90));
CHECK(geom_i.GetBeamX_pxl() == Catch::Approx(geom.GetBeamX_pxl()));
CHECK(geom_i.GetBeamY_pxl() == Catch::Approx(geom.GetBeamY_pxl()));
viewer/JFJochImageReadingWorker.cpp
+8 -5
View File
@@ -234,7 +234,7 @@ void JFJochImageReadingWorker::Analyze() {
emit imageLoaded(current_image_ptr);
}
void JFJochImageReadingWorker::FindCenter() {
void JFJochImageReadingWorker::FindCenter(const UnitCell& calibrant) {
QMutexLocker locker(&m);
if (!current_image_ptr)
return;
@@ -242,7 +242,7 @@ void JFJochImageReadingWorker::FindCenter() {
DiffractionGeometry geom = current_image_ptr->Dataset().experiment.GetDiffractionGeometry();
try {
GuessGeometry(geom, current_image_ptr->ImageData().spots, LAB6_CELL_A);
GuessGeometry(geom, current_image_ptr->ImageData().spots, calibrant);
} catch (const JFJochException &e) {
logger.ErrorException(e);
return;
@@ -258,10 +258,13 @@ void JFJochImageReadingWorker::FindCenter() {
.PoniRot2_rad(geom.GetPoniRot2_rad())
.PoniRot3_rad(geom.GetPoniRot3_rad());
UpdateDataset_i(new_experiment);
QVector<float> rings;
for (int i = 1; i < 7; i++)
rings.push_back(LAB6_CELL_A / sqrtf(i));
std::vector<float> ring_Q = CalculateXtalRings(calibrant);
QVector<float> rings;
for (int i = 0; i < 6 && i < ring_Q.size(); i++) {
rings.push_back(2 * M_PI / ring_Q[i]);
}
emit setRings(rings);
}
viewer/JFJochImageReadingWorker.h
+2 -1
View File
@@ -22,6 +22,7 @@ Q_DECLARE_METATYPE(std::shared_ptr<const JFJochReaderImage>)
Q_DECLARE_METATYPE(DiffractionExperiment)
Q_DECLARE_METATYPE(SpotFindingSettings)
Q_DECLARE_METATYPE(IndexingSettings)
Q_DECLARE_METATYPE(UnitCell)
class JFJochImageReadingWorker : public QObject {
Q_OBJECT
@@ -77,7 +78,7 @@ public slots:
void UpdateDataset(const DiffractionExperiment& experiment);
void FindCenter();
void FindCenter(const UnitCell& calibrant);
void Analyze();
void UpdateSpotFindingSettings(const SpotFindingSettings &settings, const IndexingSettings &indexing, int64_t max_spots, bool reanalyze);
viewer/JFJochViewerSidePanel.cpp
+47 -2
View File
@@ -119,8 +119,15 @@ JFJochViewerSidePanel::JFJochViewerSidePanel(QWidget *parent) : QWidget(parent)
connect(analyzeButton, &QPushButton::clicked,[this] {emit analyze();});
layout->addWidget(analyzeButton);
auto findBeamCenterButton = new QPushButton("Calibrate detector (LaB6)", this);
connect(findBeamCenterButton, &QPushButton::clicked,[this] {emit findBeamCenter();});
// Calibrant selection combo box
layout->addWidget(new TitleLabel("Powder geometry calibration", this));
calibrantCombo = new QComboBox(this);
layout->addWidget(calibrantCombo);
updateCalibrantList();
auto findBeamCenterButton = new QPushButton("Calibrate detector", this);
connect(findBeamCenterButton, &QPushButton::clicked,this, &JFJochViewerSidePanel::findBeamCenterClicked);
layout->addWidget(findBeamCenterButton);
// Add preset ice rings button below LaB6 calibration
@@ -136,6 +143,40 @@ JFJochViewerSidePanel::JFJochViewerSidePanel(QWidget *parent) : QWidget(parent)
setLayout(layout); // Set the layout to the widget
}
void JFJochViewerSidePanel::updateCalibrantList() {
calibrantCombo->clear();
calibrantCombo->addItem("LaB6");
calibrantCombo->addItem("Silver Behenate");
if (sample_cell)
calibrantCombo->addItem(QString("Current sample (%1 %2 %3 %4 %5 %6)")
.arg(QString::number(sample_cell->a, 'f', 1))
.arg(QString::number(sample_cell->b, 'f', 1))
.arg(QString::number(sample_cell->c, 'f', 1))
.arg(QString::number(sample_cell->alpha, 'f', 1))
.arg(QString::number(sample_cell->beta, 'f', 1))
.arg(QString::number(sample_cell->gamma, 'f', 1)));
calibrantCombo->setCurrentIndex(0);
}
void JFJochViewerSidePanel::findBeamCenterClicked() {
UnitCell uc(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90, 90, 90);
switch (calibrantCombo->currentIndex()) {
case 1:
// T. C. Huang , H. Toraya, T. N. Blanton, Y. Wu, J. Appl. Cryst. 26 (1993), 180-184.
uc = UnitCell(5.1769, 4.7218, 58.380, 89.440, 89.634, 75.854);
break;
case 2:
if (sample_cell)
uc = sample_cell.value();
break;
default:
break;
}
emit findBeamCenter(uc);
}
void JFJochViewerSidePanel::setRings(const QVector <float> &v) {
QString txt;
for (float i : v) {
@@ -192,6 +233,10 @@ void JFJochViewerSidePanel::editingFinished() {
}
void JFJochViewerSidePanel::loadImage(std::shared_ptr<const JFJochReaderImage> image) {
if (image)
sample_cell = image->Dataset().experiment.GetUnitCell();
updateCalibrantList();
emit imageLoaded(image);
}
viewer/JFJochViewerSidePanel.h
+7 -1
View File
@@ -16,10 +16,13 @@
class JFJochViewerSidePanel : public QWidget {
Q_OBJECT
std::optional<UnitCell> sample_cell;
QLineEdit *res_rings_edit = nullptr;
QCheckBox *autoResRingsCheckBox;
QCheckBox *resRingsCheckBox;
QComboBox* calibrantCombo{nullptr}; // stores current calibrant selection
JFJochViewerSidePanelChart *chart = nullptr;
signals:
void showSpots(bool input);
@@ -31,7 +34,7 @@ signals:
void showSaturatedPixels(bool input);
void showPredictions(bool input);
void highlightIceRings(bool input);
void findBeamCenter();
void findBeamCenter(const UnitCell &input);
void analyze();
void imageLoaded(std::shared_ptr<const JFJochReaderImage> image);
@@ -48,6 +51,9 @@ private slots:
void highlightIceRingsToggled(bool input);
void saturatedPixelsToggled(bool input);
void setRings(const QVector<float> &v);
void findBeamCenterClicked();
void updateCalibrantList();
};