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

portability: replace M_PI / std::numbers::pi with a host+device-safe PI

Browse Source 
M_PI is a POSIX <math.h> extension that MSVC does not define without
_USE_MATH_DEFINES. std::numbers::pi (introduced in the viewer guard commit) is
C++20, but CUDA here is compiled as C++17 (CMAKE_CUDA_STANDARD 17) and several
common/ headers are pulled into .cu device translation units, so std::numbers is
not available there.

Add common/JFJochMath.h with a dependency-free `constexpr double PI` that works
in host code (including MSVC), in CUDA device code, and under C++17/20, and use
it everywhere:
- common/ and image_analysis/ (incl. CUDA .cu): 78 M_PI occurrences, 22 files
- broker/OpenAPIConvert.cpp
- viewer/: the 5 files that used std::numbers::pi now use PI, for one consistent
  convention across the codebase

Verified to build: JFJochImageAnalysis (incl. CUDA), jfjoch_viewer, JFJochBroker.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-18 14:54:01 +02:00
co-authored by Claude Opus 4.8
parent d0a1175c65
commit 5d9e1be814
29 changed files with 120 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files
broker/OpenAPIConvert.cpp
+3 -2
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../common/JFJochMath.h"
#include "OpenAPIConvert.h"
// From https://en.cppreference.com/w/cpp/string/byte/tolower
@@ -460,8 +461,8 @@ ROIDefinition Convert(const org::openapitools::server::model::Roi_definitions& i
output.circles.emplace_back(ROICircle(i.getName(), i.getCenterXPxl(), i.getCenterYPxl(), i.getRadiusPxl()));
for (const auto &i: input.getAzim().getRois())
output.azimuthal.emplace_back(ROIAzimuthal(i.getName(),
(i.getQMaxRecipA() == 0.0) ? 0.0 : 2.0f * M_PI / i.getQMaxRecipA(),
(i.getQMinRecipA() == 0.0) ? 0.0 : 2.0f * M_PI / i.getQMinRecipA()));
(i.getQMaxRecipA() == 0.0) ? 0.0 : 2.0f * PI / i.getQMaxRecipA(),
(i.getQMinRecipA() == 0.0) ? 0.0 : 2.0f * PI / i.getQMinRecipA()));
return output;
}
common/AzimuthalIntegrationMapping.cpp
+6 -5
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include <cmath>
#include <thread>
#include <future>
@@ -145,8 +146,8 @@ void AzimuthalIntegrationMapping::SetupPixel(const DiffractionGeometry &geom,
corrections[pxl] = corr;
if (d > 0) {
float q = 2.0f * static_cast<float>(M_PI) / d;
pixel_to_bin[pxl] = settings.GetBin(q, phi_rad * 180.0 / M_PI);
float q = 2.0f * static_cast<float>(PI) / d;
pixel_to_bin[pxl] = settings.GetBin(q, phi_rad * 180.0 / PI);
}
}
@@ -187,9 +188,9 @@ void AzimuthalIntegrationMapping::UpdateMaxBinNumber() {
for (int j = 0; j < settings.GetAzimuthalBinCount(); j++) {
for (int i = 0; i < settings.GetQBinCount(); i++) {
bin_to_q[j * settings.GetQBinCount() + i] = static_cast<float>(settings.GetQSpacing_recipA() * (i + 0.5) + settings.GetLowQ_recipA());
bin_to_d[j * settings.GetQBinCount() + i] = 2.0f * static_cast<float>(M_PI) / bin_to_q[j * settings.GetQBinCount() + i];
bin_to_2theta[j * settings.GetQBinCount() + i] = 2.0f * asinf(bin_to_q[i] * wavelength / (4.0f * static_cast<float>(M_PI))) * 180.0f /
static_cast<float>(M_PI);
bin_to_d[j * settings.GetQBinCount() + i] = 2.0f * static_cast<float>(PI) / bin_to_q[j * settings.GetQBinCount() + i];
bin_to_2theta[j * settings.GetQBinCount() + i] = 2.0f * asinf(bin_to_q[i] * wavelength / (4.0f * static_cast<float>(PI))) * 180.0f /
static_cast<float>(PI);
bin_to_phi[j * settings.GetQBinCount() + i] = static_cast<float>(j) * 360.0f / static_cast<float>(settings.GetAzimuthalBinCount());
}
}
common/AzimuthalIntegrationSettings.h
+3 -2
View File
@@ -3,6 +3,7 @@
#pragma once
#include "JFJochMath.h"
#include <optional>
#include <cstdint>
@@ -11,8 +12,8 @@ class AzimuthalIntegrationSettings {
bool polarization_correction = true;
float high_q_recipA = 5.0;
float low_q_recipA = 0.1;
float bkg_estimate_high_q_recipA = 2.0f * M_PI / 3.0;
float bkg_estimate_low_q_recipA = 2.0f * M_PI / 5.0;
float bkg_estimate_high_q_recipA = 2.0f * PI / 3.0;
float bkg_estimate_low_q_recipA = 2.0f * PI / 5.0;
float q_spacing = 0.05;
int32_t azim_bins = 1;
common/Coord.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include <cmath>
#include "Coord.h"
#include "JFJochException.h"
@@ -143,7 +144,7 @@ float angle_deg(const Coord &c1, const Coord &c2) {
else if (cos_ang < -1.0f && cos_ang > -1.0f - epsilon)
cos_ang = -1.0f;
return acosf(cos_ang) * (180.0f / static_cast<float>(M_PI));
return acosf(cos_ang) * (180.0f / static_cast<float>(PI));
}
void Coord::swap(Coord &other) noexcept {
common/CrystalLattice.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include <cmath>
#include "CrystalLattice.h"
#include "JFJochException.h"
@@ -9,7 +10,7 @@
#include "gemmi/unitcell.hpp"
#include "gemmi/cellred.hpp"
#define DEG_TO_RAD static_cast<float>(M_PI/180.0)
#define DEG_TO_RAD static_cast<float>(PI/180.0)
CrystalLattice::CrystalLattice(const UnitCell &cell) {
vec[0] = {cell.a, 0, 0};
common/DiffractionGeometry.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include <cmath>
#include "DiffractionGeometry.h"
@@ -48,7 +49,7 @@ float DiffractionGeometry::Phi_rad(float x, float y) const {
auto lab = LabCoord(x, y);
auto v = atan2f(lab.y, lab.x);
if (v < 0)
v += 2.0f * M_PI;
v += 2.0f * PI;
return v;
}
@@ -58,7 +59,7 @@ float DiffractionGeometry::PxlToRes(float x, float y) const {
}
float DiffractionGeometry::PxlToQ(float x, float y) const {
return 2.0f * M_PI / PxlToRes(x,y);
return 2.0f * PI / PxlToRes(x,y);
}
float DiffractionGeometry::PxlToRes(float dist_pxl) const {
@@ -87,7 +88,7 @@ float DiffractionGeometry::DistFromEwaldSphere(const Coord &recip) const {
}
float DiffractionGeometry::CalcAzIntSolidAngleCorr(float q) const {
float sin_theta = q * wavelength_A / (4 * static_cast<float>(M_PI));
float sin_theta = q * wavelength_A / (4 * static_cast<float>(PI));
float cos_2theta = 1.0f - 2.0f * sin_theta * sin_theta; // cos(2*alpha) = 1 - 2 * sin(alpha)^2
float cos_2theta_3 = cos_2theta * cos_2theta * cos_2theta;
return cos_2theta_3;
common/GoniometerAxis.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include <cmath>
#include <algorithm>
@@ -138,6 +139,6 @@ std::vector<double> GoniometerAxis::GetAngleContainerEnd(int64_t max_image_numbe
}
RotMatrix GoniometerAxis::GetTransformationAngle(float angle_deg) const {
auto angle_rad = angle_deg / 180.0f * static_cast<float>(M_PI);
auto angle_rad = angle_deg / 180.0f * static_cast<float>(PI);
return {angle_rad, axis};
}
common/JFJochMath.h
+10
View File
@@ -0,0 +1,10 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
// pi as a plain constexpr so it is usable everywhere:
// - host code on MSVC, which does not define M_PI without _USE_MATH_DEFINES;
// - CUDA device code, which is compiled as C++17 here, so PI
// (C++20) is not available.
constexpr double PI = 3.14159265358979323846;
common/ROIAzimuthal.cpp
+3 -2
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "JFJochMath.h"
#include "ROIAzimuthal.h"
#include "JFJochException.h"
@@ -31,11 +32,11 @@ bool ROIAzimuthal::CheckROI(int64_t x, int64_t y, float resolution) const {
}
float ROIAzimuthal::GetQMax_recipA() const {
return 2.0f * M_PI / d_min_A;
return 2.0f * PI / d_min_A;
}
float ROIAzimuthal::GetQMin_recipA() const {
return 2.0f * M_PI / d_max_A;
return 2.0f * PI / d_max_A;
}
ROIConfig ROIAzimuthal::ExportMetadata() const {
image_analysis/bragg_prediction/BraggPrediction.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "BraggPrediction.h"
#include "../bragg_integration/SystematicAbsence.h"
@@ -46,7 +47,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
const float epsilon = 1e-5f;
const float s0_sq = S0 * S0;
const float rad_to_deg = 180.0f / static_cast<float>(M_PI);
const float rad_to_deg = 180.0f / static_cast<float>(PI);
int i = 0;
image_analysis/bragg_prediction/BraggPredictionGPU.cu
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "BraggPredictionGPU.h"
#ifdef JFJOCH_USE_CUDA
@@ -16,7 +17,7 @@ namespace {
__device__ inline float angle_from_ewald_sphere_deg(const Coord &S0, float recip_x, float recip_y, float recip_z, float recip_sq) {
const float epsilon = 1e-5f;
const float rad_to_deg = 180.0f / static_cast<float>(M_PI);
const float rad_to_deg = 180.0f / static_cast<float>(PI);
const float s0_sq = S0.x * S0.x + S0.y * S0.y + S0.z * S0.z;
const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z;
image_analysis/bragg_prediction/BraggPredictionRot.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "BraggPredictionRot.h"
#include "../bragg_integration/SystematicAbsence.h"
@@ -46,8 +47,8 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
int i = 0;
const float mos_angle_rad = settings.mosaicity_deg * static_cast<float>(M_PI) / 180.f;
const float half_wedge_angle_rad = settings.wedge_deg * static_cast<float>(M_PI) / 180.f / 2.0f ;
const float mos_angle_rad = settings.mosaicity_deg * static_cast<float>(PI) / 180.f;
const float half_wedge_angle_rad = settings.wedge_deg * static_cast<float>(PI) / 180.f / 2.0f ;
for (int h = -settings.max_hkl; h <= settings.max_hkl; h++) {
// Precompute A* h contribution
@@ -134,7 +135,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
.h = h,
.k = k,
.l = l,
.delta_phi_deg = phi * 180.0f / static_cast<float>(M_PI),
.delta_phi_deg = phi * 180.0f / static_cast<float>(PI),
.predicted_x = x,
.predicted_y = y,
.observed_x = NAN,
image_analysis/bragg_prediction/BraggPredictionRotGPU.cu
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "BraggPredictionRotGPU.h"
#ifdef JFJOCH_USE_CUDA
@@ -150,7 +151,7 @@ namespace {
out[count].h = h;
out[count].k = k;
out[count].l = l;
out[count].delta_phi_deg = phi * 180.0 / M_PI;
out[count].delta_phi_deg = phi * 180.0 / PI;
out[count].predicted_x = x;
out[count].predicted_y = y;
out[count].observed_x = NAN;
@@ -209,8 +210,8 @@ namespace {
kc.one_over_wavelength = 1.0f / geom.GetWavelength_A();
// Store mosaicity and wedge in radians for partiality calculation
kc.mos_angle_rad = settings.mosaicity_deg * static_cast<float>(M_PI) / 180.0f;
kc.wedge_angle_rad = settings.wedge_deg * static_cast<float>(M_PI) / 180.0f;
kc.mos_angle_rad = settings.mosaicity_deg * static_cast<float>(PI) / 180.0f;
kc.wedge_angle_rad = settings.wedge_deg * static_cast<float>(PI) / 180.0f;
kc.min_zeta = settings.min_zeta;
kc.mosaicity_multiplier = settings.mosaicity_multiplier;
image_analysis/geom_refinement/AssignSpotsToRings.cpp
+5 -4
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "AssignSpotsToRings.h"
#include <vector>
@@ -163,7 +164,7 @@ std::vector<float> CalculateXtalRings(const UnitCell &cell, int hkl_max) {
for (int l = 0; l <= hkl_max; l++) {
if (h == 0 && k == 0 && l == 0) continue;
auto p = Astar * h + Bstar * k + Cstar * l;
float Q = 2.0f * M_PI * p.Length();
float Q = 2.0f * PI * p.Length();
u.push_back(Q);
}
}
@@ -222,7 +223,7 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v,
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]);
auto cluster_annot = GuessInitialGeometry(geom, v, 2 * PI / ring_q[0]);
std::vector<RingOptimizerInput> optimizer_input;
@@ -231,7 +232,7 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v,
while (cluster_idx < cluster_annot.size() && ring_idx < ring_q.size()) {
float obs_q = 2 * M_PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs);
float obs_q = 2 * PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs);
if (std::fabs(ring_q[ring_idx] - obs_q) < 0.1) {
std::cout << "Found ring " << ring_idx
<< " with q_ring = " << ring_q[ring_idx]
@@ -263,7 +264,7 @@ void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &
std::vector<RingOptimizerInput> optimizer_input;
for (const auto& s: v) {
float q_obs = 2 * M_PI / geom.PxlToRes(s.x, s.y);
float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y);
for (const auto &q : q_ring) {
if (std::fabs(q - q_obs) < 0.1) {
optimizer_input.push_back({s.x, s.y, q});
image_analysis/geom_refinement/LatticeReduction.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "LatticeReduction.h"
#include <Eigen/Dense>
@@ -251,7 +252,7 @@ CrystalLattice AngleAxisAndCellToLattice(const double rod[3],
CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex) {
return AngleAxisAndCellToLattice(rod, lengths,
/*alpha=*/M_PI / 2.0,
/*beta =*/M_PI / 2.0,
/*gamma=*/hex ? 2.0 * M_PI / 3.0 : M_PI / 2.0);
/*alpha=*/PI / 2.0,
/*beta =*/PI / 2.0,
/*gamma=*/hex ? 2.0 * PI / 3.0 : PI / 2.0);
}
image_analysis/geom_refinement/RingOptimizer.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "RingOptimizer.h"
#include "ceres/ceres.h"
@@ -11,7 +12,7 @@ struct RingResidual {
: obs_x(x), obs_y(y),
lambda(lambda),
pixel_size(pixel_size),
expected_len_recip_sq(expected_q * expected_q / (4.0 * M_PI * M_PI)) {}
expected_len_recip_sq(expected_q * expected_q / (4.0 * PI * PI)) {}
template<typename T>
bool operator()(const T* const center_x, const T* const center_y,
image_analysis/geom_refinement/XtalOptimizer.cpp
+15 -14
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include <Eigen/Dense>
#include "XtalOptimizer.h"
@@ -293,9 +294,9 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
// Triclinic: initialize a,b,c and α,β,γ from current unit cell
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
auto uc = data.latt.GetUnitCell();
latt_vec2[0] = uc.alpha * M_PI / 180.0;
latt_vec2[1] = uc.beta * M_PI / 180.0;
latt_vec2[2] = uc.gamma * M_PI / 180.0;
latt_vec2[0] = uc.alpha * PI / 180.0;
latt_vec2[1] = uc.beta * PI / 180.0;
latt_vec2[2] = uc.gamma * PI / 180.0;
break;
}
@@ -316,7 +317,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
if (data.axis) {
const float angle_deg = data.axis->GetAngle_deg(i) + data.axis->GetWedge_deg() / 2.0;
angle_rad = angle_deg * M_PI / 180.0;
angle_rad = angle_deg * PI / 180.0;
rot_matr = data.axis->GetTransformationAngle(angle_deg);
}
@@ -380,7 +381,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
if (!data.refine_detector_angles) {
problem.SetParameterBlockConstant(detector_rot);
} else {
const double rot_range = 3.0 / 180.0 * M_PI;
const double rot_range = 3.0 / 180.0 * PI;
for (int i = 0; i < 2; ++i) {
problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rot_range);
problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rot_range);
@@ -403,14 +404,14 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
}
if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) {
const double beta_lo = std::max(1e-6, M_PI * (data.min_angle_deg / 180.0));
const double beta_hi = std::min(M_PI - 1e-6, M_PI * (data.max_angle_deg / 180.0));
const double beta_lo = std::max(1e-6, PI * (data.min_angle_deg / 180.0));
const double beta_hi = std::min(PI - 1e-6, PI * (data.max_angle_deg / 180.0));
problem.SetParameterLowerBound(latt_vec2, 0, beta_lo);
problem.SetParameterUpperBound(latt_vec2, 0, beta_hi);
} else if (data.crystal_system == gemmi::CrystalSystem::Triclinic) {
// α, β, γ bounds (radians)
const double alo = M_PI * (data.min_angle_deg / 180.0);
const double ahi = M_PI * (data.max_angle_deg / 180.0);
const double alo = PI * (data.min_angle_deg / 180.0);
const double ahi = PI * (data.max_angle_deg / 180.0);
for (int i = 0; i < 3; ++i) {
problem.SetParameterLowerBound(latt_vec2, i, alo);
problem.SetParameterUpperBound(latt_vec2, i, ahi);
@@ -450,19 +451,19 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
data.axis.value().Axis(Coord(rot_vec[0], rot_vec[1], rot_vec[2]));
if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic)
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_PI / 2.0);
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) {
latt_vec1[1] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_PI / 2.0);
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
} else if (data.crystal_system == gemmi::CrystalSystem::Cubic) {
latt_vec1[1] = latt_vec1[0];
latt_vec1[2] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_PI / 2.0);
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
} else if (data.crystal_system == gemmi::CrystalSystem::Hexagonal) {
latt_vec1[1] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0);
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,PI / 2.0, PI / 2.0, 2.0 * PI / 3.0);
} else if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) {
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, latt_vec2[0], M_PI / 2.0);
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, latt_vec2[0], PI / 2.0);
} else {
// Triclinic via the same generic builder
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, latt_vec2[0], latt_vec2[1], latt_vec2[2]);
image_analysis/indexing/AnalyzeIndexing.cpp
+6 -5
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include <cstdint>
#include <vector>
#include <cmath>
@@ -20,11 +21,11 @@ namespace {
}
inline float deg_to_rad(float deg) {
return deg * (static_cast<float>(M_PI) / 180.0f);
return deg * (static_cast<float>(PI) / 180.0f);
}
inline float rad_to_deg(float rad) {
return rad * (180.0f / static_cast<float>(M_PI));
return rad * (180.0f / static_cast<float>(PI));
}
// Wrap to [-180, 180] (useful for residuals)
@@ -139,7 +140,7 @@ namespace {
float s1 = phi_ref + delta;
float s2 = phi_ref - delta;
const float two_pi = 2.0f * static_cast<float>(M_PI);
const float two_pi = 2.0f * static_cast<float>(PI);
auto shift_near = [&](float x) {
if (!std::isfinite(x))
return std::numeric_limits<float>::quiet_NaN();
@@ -153,9 +154,9 @@ namespace {
return std::numeric_limits<float>::quiet_NaN();
// fmod can return negative values; normalize to [-π, π]
if (shifted < -static_cast<float>(M_PI))
if (shifted < -static_cast<float>(PI))
shifted += two_pi;
else if (shifted > static_cast<float>(M_PI))
else if (shifted > static_cast<float>(PI))
shifted -= two_pi;
return shifted + span_center;
image_analysis/indexing/FFTIndexer.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "FFTIndexer.h"
#include <Eigen/Eigen>
#include "PostIndexingRefinement.h"
@@ -23,7 +24,7 @@ FFTIndexer::FFTIndexer(const IndexingSettings &settings)
// from when the extent was mistakenly written as 2*pi/d (the Q convention); it is kept
// because the exact amount sets which marginal frames index - rounding it to a nearby
// integer shifts the indexing rate ~0.5-1% (validated on the lyso datasets).
const float oversampling = 2.0f * static_cast<float>(M_PI);
const float oversampling = 2.0f * static_cast<float>(PI);
const float one_over_d_max = 1.0f / settings.GetFFT_HighResolution_A();
histogram_spacing = 1.0f / (2.0f * max_length_A);
@@ -52,7 +53,7 @@ void FFTIndexer::SetupDirectionVectors() {
direction_vectors.reserve(static_cast<size_t>(nDirections));
const double phi = (1.0 + std::sqrt(5.0)) / 2.0; // Golden ratio
const double golden_angle = 2.0 * M_PI / phi;
const double golden_angle = 2.0 * PI / phi;
for (int i = 0; i < nDirections; i++) {
// Half-sphere distribution (z in [0,1])
@@ -148,7 +149,7 @@ std::vector<Coord> FFTIndexer::FilterFFTResults() const {
std::vector<Coord> ret;
// Remove vectors less than 5 deg apart, as most likely these are colinear
const float COS_5_DEG = std::cos(5.0f * M_PI / 180.0f);
const float COS_5_DEG = std::cos(5.0f * PI / 180.0f);
// Minimum relative amplitude to accept a shorter vector (fundamental)
// over a longer one (harmonic).
image_analysis/indexing/PostIndexingRefinement.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "PostIndexingRefinement.h"
#include <iostream>
@@ -175,9 +176,9 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
}
// Filter for wrong angles
float alpha = std::acos(cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
float beta = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / M_PI;
float gamma = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized())) * 180.0f / M_PI;
float alpha = std::acos(cell_rows.row(1).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / PI;
float beta = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(2).normalized())) * 180.0f / PI;
float gamma = std::acos(cell_rows.row(0).normalized().dot(cell_rows.row(1).normalized())) * 180.0f / PI;
if (alpha < p.min_angle_deg || alpha > p.max_angle_deg ||
beta < p.min_angle_deg || beta > p.max_angle_deg ||
image_analysis/lattice_search/LatticeSearch.cpp
+4 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "LatticeSearch.h"
#include <gemmi/cellred.hpp>
@@ -315,9 +316,9 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
if (c.cond_BC && fabs((uc_reduced.b - uc_reduced.c) / (0.5 * (uc_reduced.b + uc_reduced.c))) > dist_tolerance)
ok = false;
double expected_alpha = acos(c.cond_D / sqrt(B*C)) * 180 / M_PI;
double expected_beta = acos(c.cond_E / sqrt(A*C)) * 180 / M_PI;
double expected_gamma = acos(c.cond_F / sqrt(A*B)) * 180 / M_PI;
double expected_alpha = acos(c.cond_D / sqrt(B*C)) * 180 / PI;
double expected_beta = acos(c.cond_E / sqrt(A*C)) * 180 / PI;
double expected_gamma = acos(c.cond_F / sqrt(A*B)) * 180 / PI;
if (fabs(expected_alpha - uc_reduced.alpha) > angle_tolerance)
ok = false;
image_analysis/pixel_refinement/PixelRefine.cpp
+3 -2
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "PixelRefine.h"
#include <algorithm>
@@ -499,7 +500,7 @@ void PixelRefine::Run(const T *image,
beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2,
q_sq, eps_r, eps_t_sq))
continue;
const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (M_PI * R1 * R1);
const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (PI * R1 * R1);
const double R0_eff_sq = R0 * R0 + g.R_bw_sq;
const double P_rad = std::exp(-eps_r * eps_r / R0_eff_sq);
const double v = std::max(g.Ibkg, 1.0);
@@ -597,7 +598,7 @@ void PixelRefine::Run(const T *image,
// Tangential profile shape (area-normalized) -> the fit template, using the
// per-reflection R1_eff (Term 2).
const double R1 = g.R1_eff;
const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (M_PI * R1 * R1);
const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (PI * R1 * R1);
prof_full += P_t; // whole shoebox, on- or off-detector
// Only real, unmasked detector pixels carry signal.
image_analysis/rotation_indexer/RotationIndexer.cpp
+2 -1
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "RotationIndexer.h"
#include "../geom_refinement/XtalOptimizer.h"
@@ -84,7 +85,7 @@ void RotationIndexer::RunIndexing() {
break;
// Ignore lattices oriented by less than 3.0 degree
if (l.rotation_vector.Length() < 3.0 * M_PI / 180.0)
if (l.rotation_vector.Length() < 3.0 * PI / 180.0)
continue;
RotMatrix rot(l.rotation_vector.Length(), l.rotation_vector.Normalize());
image_analysis/spot_finding/SpotUtils.cpp
+5 -4
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "../../common/JFJochMath.h"
#include "SpotUtils.h"
#include "../../common/ResolutionShells.h"
@@ -25,10 +26,10 @@ void MarkIceRings(std::vector<SpotToSave> &spots, float tolerance_q_recipA) {
std::vector<float> ice_rings_q;
for (const auto &i: ICE_RING_RES_A)
ice_rings_q.push_back(2 * M_PI / i);
ice_rings_q.push_back(2 * PI / i);
for (auto &s: spots) {
auto spot_q = 2 * M_PI / s.d_A;
auto spot_q = 2 * PI / s.d_A;
bool tmp = false;
for (const auto &q: ice_rings_q)
tmp |= (fabs(spot_q - q) < tolerance_q_recipA);
@@ -65,8 +66,8 @@ void FilterSpuriousHighResolutionSpots(std::vector<SpotToSave> &spots, float thr
// Avoid division by zero; d_A should be > 0 in valid data
if (d1 <= 0.0f || d2 <= 0.0f)
continue;
float q1 = 2 * M_PI / d1;
float q2 = 2 * M_PI / d2;
float q1 = 2 * PI / d1;
float q2 = 2 * PI / d2;
float dq = q2 - q1; // should be >= 0 due to sorting
if (dq > threshold) {
cut_index = i + 1; // keep up to i inclusive
viewer/JFJochImageReadingWorker.cpp
+2 -2
View File
@@ -3,7 +3,7 @@
#include <cerrno>
#include <cstring>
#include <numbers>
#include "../common/JFJochMath.h"
#ifndef _WIN32
#include <sys/stat.h>
#include <fcntl.h>
@@ -539,7 +539,7 @@ void JFJochImageReadingWorker::FindCenter(const UnitCell& calibrant, bool guess)
QVector<float> rings;
for (int i = 0; i < 15 && i < ring_Q.size(); i++) {
rings.push_back(2 * std::numbers::pi / ring_Q[i]);
rings.push_back(2 * PI / ring_Q[i]);
}
emit setRings(rings);
}
viewer/image_viewer/JFJochAzIntImage.cpp
+2 -2
View File
@@ -5,7 +5,7 @@
#include <QVBoxLayout>
#include "../../common/JFJochException.h"
#include <numbers>
#include "../../common/JFJochMath.h"
JFJochAzIntImage::JFJochAzIntImage(QWidget *parent)
: JFJochImage(parent) {
@@ -110,7 +110,7 @@ void JFJochAzIntImage::mouseDoubleClickEvent(QMouseEvent *event) {
float phi = image->Dataset().az_int_bin_to_phi[idx];
auto geom = image->Dataset().experiment.GetDiffractionGeometry();
auto coord = geom.ResPhiToPxl(2 * std::numbers::pi / q, phi / 180.0 * std::numbers::pi);
auto coord = geom.ResPhiToPxl(2 * PI / q, phi / 180.0 * PI);
emit zoomOnBin(QPointF(coord.first, coord.second));
}
}
viewer/image_viewer/JFJochDiffractionImage.cpp
+8 -8
View File
@@ -3,7 +3,7 @@
#include "JFJochDiffractionImage.h"
#include "../../common/DiffractionGeometry.h"
#include <numbers>
#include "../../common/JFJochMath.h"
#include <QGraphicsPixmapItem>
#include <QGraphicsSimpleTextItem>
@@ -240,9 +240,9 @@ void JFJochDiffractionImage::DrawResolutionRings() {
continue;
auto [x1,y1] = geom.ResPhiToPxl(d, 0);
auto [x2,y2] = geom.ResPhiToPxl(d, std::numbers::pi / 2);
auto [x3,y3] = geom.ResPhiToPxl(d, std::numbers::pi);
auto [x4,y4] = geom.ResPhiToPxl(d, 3.0 * std::numbers::pi / 2);
auto [x2,y2] = geom.ResPhiToPxl(d, PI / 2);
auto [x3,y3] = geom.ResPhiToPxl(d, PI);
auto [x4,y4] = geom.ResPhiToPxl(d, 3.0 * PI / 2);
auto x_min = std::min({x1, x2, x3, x4});
auto x_max = std::max({x1, x2, x3, x4});
@@ -253,9 +253,9 @@ void JFJochDiffractionImage::DrawResolutionRings() {
addOverlayItem(scene()->addEllipse(boundingRect, pen));
auto [x5,y5] = geom.ResPhiToPxl(d, phi_offset + 0);
auto [x6,y6] = geom.ResPhiToPxl(d, phi_offset + std::numbers::pi / 2);
auto [x7,y7] = geom.ResPhiToPxl(d, phi_offset + std::numbers::pi);
auto [x8,y8] = geom.ResPhiToPxl(d, phi_offset + 3.0 * std::numbers::pi / 2);
auto [x6,y6] = geom.ResPhiToPxl(d, phi_offset + PI / 2);
auto [x7,y7] = geom.ResPhiToPxl(d, phi_offset + PI);
auto [x8,y8] = geom.ResPhiToPxl(d, phi_offset + 3.0 * PI / 2);
QPointF point_1(x5, y5);
QPointF point_2(x6, y6);
@@ -286,7 +286,7 @@ void JFJochDiffractionImage::DrawResolutionRings() {
scene()->addItem(textItem);
addOverlayItem(textItem);
}
phi_offset += 4.0 / 180.0 * std::numbers::pi;
phi_offset += 4.0 / 180.0 * PI;
}
}
viewer/widgets/PowderCalibrationWidget.cpp
+2 -2
View File
@@ -5,7 +5,7 @@
#include <QPushButton>
#include <QGridLayout>
#include <QLabel>
#include <numbers>
#include "../../common/JFJochMath.h"
#include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
@@ -26,7 +26,7 @@ PowderCalibrationWidget::PowderCalibrationWidget(QWidget *parent) : QWidget(pare
std::vector<float> rings = CalculateXtalRings(GetCalibrant(), 10);
QVector<float> q_rings;
for (float ring : rings) {
q_rings.append(2 * std::numbers::pi / ring);
q_rings.append(2 * PI / ring);
}
emit ringsFromCalibration(q_rings);
});
viewer/windows/JFJochViewerMetadataWindow.cpp
+5 -5
View File
@@ -10,7 +10,7 @@
#include <QMessageBox>
#include "JFJochViewerMetadataWindow.h"
#include <numbers>
#include "../../common/JFJochMath.h"
JFJochViewerMetadataWindow::JFJochViewerMetadataWindow(QWidget *parent) : JFJochHelperWindow(parent) {
@@ -136,8 +136,8 @@ void JFJochViewerMetadataWindow::datasetUpdate() {
tmp.BeamX_pxl(beam_center_x->value());
tmp.BeamY_pxl(beam_center_y->value());
tmp.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A->value());
tmp.PoniRot1_rad(poni_rot1_deg->value() * std::numbers::pi / 180.0);
tmp.PoniRot2_rad(poni_rot2_deg->value() * std::numbers::pi / 180.0);
tmp.PoniRot1_rad(poni_rot1_deg->value() * PI / 180.0);
tmp.PoniRot2_rad(poni_rot2_deg->value() * PI / 180.0);
tmp.DetectIceRings(detect_ice_rings->isChecked());
if (unit_cell_enabled->isChecked()) {
@@ -177,8 +177,8 @@ void JFJochViewerMetadataWindow::datasetLoaded(std::shared_ptr<const JFJochReade
beam_center_x->setValue(dataset->experiment.GetBeamX_pxl());
beam_center_y->setValue(dataset->experiment.GetBeamY_pxl());
wavelength_A->setValue(dataset->experiment.GetWavelength_A());
poni_rot1_deg->setValue(dataset->experiment.GetPoniRot1_rad() * 180.0 / std::numbers::pi);
poni_rot2_deg->setValue(dataset->experiment.GetPoniRot2_rad() * 180.0 / std::numbers::pi);
poni_rot1_deg->setValue(dataset->experiment.GetPoniRot1_rad() * 180.0 / PI);
poni_rot2_deg->setValue(dataset->experiment.GetPoniRot2_rad() * 180.0 / PI);
detect_ice_rings->setChecked(dataset->experiment.IsDetectIceRings());
auto unit_cell = dataset->experiment.GetUnitCell();