diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index 9babc079..0020f87b 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -9,6 +9,8 @@ #include #include +#include "../../common/CrystalLattice.h" + FindCircleCenterResult FindCircleCenter(const std::vector &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 AnalyzeClusters(const std::vector& 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 CalculateCubicXtalRings(float a, int hkl_max) { +std::vector CalculateXtalRings(const UnitCell &cell, int hkl_max) { + CrystalLattice latt(cell); + + Coord Astar = latt.Astar(); + Coord Bstar = latt.Bstar(); + Coord Cstar = latt.Cstar(); + std::vector 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 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 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 GuessInitialGeometry(DiffractionGeometry &geom, const return cluster_annot; } -void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A) { - auto cluster_annot = GuessInitialGeometry(geom, v, calibrant_a_A); +void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant_uc) { + std::vector 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 ring_q = CalculateCubicXtalRings(calibrant_a_A); std::vector optimizer_input; int ring_idx = 0; @@ -246,8 +257,8 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, geom = optimizer.Run(optimizer_input); } -void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A) { - std::vector q_ring = CalculateCubicXtalRings(calibrant_a_A); +void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant) { + std::vector q_ring = CalculateXtalRings(calibrant); std::vector optimizer_input; diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index c833bf5c..835f9e59 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -5,7 +5,7 @@ #define JFJOCH_ASSIGNSPOTSTORINGS_H #include -#include +#include "../../common/UnitCell.h" #include "../../common/SpotToSave.h" #include "RingOptimizer.h" @@ -29,11 +29,12 @@ FindCircleCenterResult FindCircleCenter(const std::vector &v, std::vector> ClusterSpotsIntoRings(const std::vector& r, float eps = 1.5, int minPts = 5); std::vector AnalyzeClusters(const std::vector& r, const std::vector> &clusters); +std::vector CalculateXtalRings(const UnitCell &cell, int hkl_max = 6); std::vector CalculateCubicXtalRings( float a, int hkl_max = 4); float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A); std::vector GuessInitialGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A); -void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A); -void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A); +void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant); +void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant); #endif //JFJOCH_ASSIGNSPOTSTORINGS_H diff --git a/tests/DetGeomCalibTest.cpp b/tests/DetGeomCalibTest.cpp index 82c79c39..7328bf87 100644 --- a/tests/DetGeomCalibTest.cpp +++ b/tests/DetGeomCalibTest.cpp @@ -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 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())); diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 58c55b89..83ba0d7c 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -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 rings; - for (int i = 1; i < 7; i++) - rings.push_back(LAB6_CELL_A / sqrtf(i)); + std::vector ring_Q = CalculateXtalRings(calibrant); + + QVector rings; + for (int i = 0; i < 6 && i < ring_Q.size(); i++) { + rings.push_back(2 * M_PI / ring_Q[i]); + } emit setRings(rings); } diff --git a/viewer/JFJochImageReadingWorker.h b/viewer/JFJochImageReadingWorker.h index a3ea25de..19269499 100644 --- a/viewer/JFJochImageReadingWorker.h +++ b/viewer/JFJochImageReadingWorker.h @@ -22,6 +22,7 @@ Q_DECLARE_METATYPE(std::shared_ptr) 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); diff --git a/viewer/JFJochViewerSidePanel.cpp b/viewer/JFJochViewerSidePanel.cpp index 221f6f9d..c87bf27f 100644 --- a/viewer/JFJochViewerSidePanel.cpp +++ b/viewer/JFJochViewerSidePanel.cpp @@ -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 &v) { QString txt; for (float i : v) { @@ -192,6 +233,10 @@ void JFJochViewerSidePanel::editingFinished() { } void JFJochViewerSidePanel::loadImage(std::shared_ptr image) { + if (image) + sample_cell = image->Dataset().experiment.GetUnitCell(); + updateCalibrantList(); + emit imageLoaded(image); } diff --git a/viewer/JFJochViewerSidePanel.h b/viewer/JFJochViewerSidePanel.h index 84bae130..2c7f8f51 100644 --- a/viewer/JFJochViewerSidePanel.h +++ b/viewer/JFJochViewerSidePanel.h @@ -16,10 +16,13 @@ class JFJochViewerSidePanel : public QWidget { Q_OBJECT + std::optional 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 image); @@ -48,6 +51,9 @@ private slots: void highlightIceRingsToggled(bool input); void saturatedPixelsToggled(bool input); void setRings(const QVector &v); + + void findBeamCenterClicked(); + void updateCalibrantList(); };