diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 31dcd9f6..fd1dec6b 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -20,6 +20,8 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: **`--mode` replaces `--azint-only` and `--scale`** (see below). * jfjoch_viewer: **"Analyze dataset" can run the detector calibration** over the whole dataset. The settings panel's MX/AzInt toggle gains a third *Calib* page carrying the calibrant (now including `ice`, matching the CLI) and the rings/spots method; the run writes a `.poni` next to the other outputs and reports the fitted PONI, tilts and distance against the header values in their own window. * Powder calibration: the ring geometry fit can now read its rings off an **azimuthally-binned profile summed over a run** instead of a spot list from one image (`RingsFromAzimuthalProfile`). A powder ring is an arc, not a set of spots, and its roundness fixes the beam centre without reference to the calibrant's d-spacings or the detector distance - the one parameter Bragg data constrain worst. The extraction window is capped at half the gap to the neighbouring ring (hexagonal ice has three rings within 0.06 1/A, which a fixed window merged into one), and where only one ring is in reach the tilts are held rather than fitted, since a single ring cannot separate them from the centre. +* Powder calibration: `.poni` files write `Poni1`/`Poni2` as (beam + 0.5) x pixel size, matching pyFAI's edge-of-sensor origin; the pixel-coordinate convention is documented in `docs/DETECTOR_GEOMETRY.md`. +* jfjoch_viewer: the beam-centre cross is aligned with the spot and prediction markers. * rugnux: De-novo **space-group search** substantially more robust - centering ranked by net absences and judged on absent-class strength, merohedral-twin over-promotion vetoed, and genuine high-symmetry groups recovered on weak data. * rugnux: The space-group search takes systematic absences from the merge of all observations, needs at least three control reflections on an axial row to claim a **screw axis**, and no longer alters the production merge. * rugnux: The space-group search no longer **starves on a low-ISa merge**: its fixed `I/sigma >= 3` cut could select nothing at all, since a merged sigma is floored so that no reflection reads above `ISa = 1/b`, leaving every operator correlation undefined and the point group at 1. The cut is now capped at the merge's own quantile, and is unchanged on healthy merges. diff --git a/docs/DETECTOR_GEOMETRY.md b/docs/DETECTOR_GEOMETRY.md index 15534556..aec81e9e 100644 --- a/docs/DETECTOR_GEOMETRY.md +++ b/docs/DETECTOR_GEOMETRY.md @@ -12,6 +12,33 @@ where beam from the sample is perpendicular to the detector surface and not to t is defined with three rotation angles: `rot1` (rotating detector right), `rot2` (rotating detector downwards), `rot3` (rotating detector clockwise). See [PyFAI documentation](https://pyfai.readthedocs.io/en/stable/) for more details. +## What a pixel coordinate means: (0, 0) is the centre of the first pixel + +Pixel coordinates in Jungfraujoch and rugnux are **0-based and pixel-centred**: an integer coordinate +is the *centre* of that pixel, so pixel *i* covers [*i* − 0.5, *i* + 0.5) and the sensor spans +−0.5 … width − 0.5. A beam centre of 948.0 × 546.0 sits in the middle of pixel [546][948], not on any +of its corners; 948.5 is the boundary between pixel 948 and 949. + +This holds throughout the code: spot and reflection centroids are intensity-weighted sums of the +integer pixel indices, the resolution and azimuthal-bin maps evaluate pixel (col, row) at exactly +(col, row), and a fractional coordinate is turned back into a pixel index by rounding, not by +truncation. The same convention applies to every coordinate the system exposes — the beam centre +(`beam_x_pxl`/`beam_y_pxl` in the API and broker configuration, `--beam-x`/`--beam-y` in rugnux, +`beam_center_x`/`beam_center_y` in NXmx and in the CBOR stream), the spot and predicted-reflection +positions written to HDF5, and the PONI reported by `--mode calibration`. + +Other programs place the origin differently, and the difference is worth half a pixel — enough to +matter when a geometry is copied between programs and then refined: + +| Convention | Beam centre equivalent to our *x* = 948.0 | +|---|---| +| Jungfraujoch, rugnux | 948.0 | +| XDS (`ORGX`/`ORGY`) | 949.0 — also pixel-centred, but pixels are numbered from 1 | +| Measured from the edge of the sensor, in length units — pyFAI (`Poni1`/`Poni2`), DIALS/dxtbx | (948.0 + 0.5) × pixel size, because the centre of pixel *i* is at (*i* + 0.5) × pixel size from the edge | + +The `.poni` file written by `rugnux --mode calibration` is in pyFAI's frame and so already carries +that half pixel; the pixel values the same run reports are ours. + ## Macromolecular crystallography convention for the vertical direction One place of confusion is the convention to have point (0,0) of the detector in the top left corner of the detector, with Y values increasing downwards. This is also consistent with computer image formats. diff --git a/docs/HDF5.md b/docs/HDF5.md index 9aaba923..a0ac8adc 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -151,7 +151,7 @@ File-level HDF5 attributes `file_name`, `file_time`, `HDF5_Version` are also set | Field | Std | Units | |-------|:---:|-------| | `depends_on` | NXmx | → `transformations/rot3` | -| `beam_center_x`, `beam_center_y` | NXmx | pixel | +| `beam_center_x`, `beam_center_y` | NXmx | pixel (0.0 = centre of the first pixel, see [DETECTOR_GEOMETRY](DETECTOR_GEOMETRY.md)) | | `distance` | NXmx | m | | `count_time`, `frame_time` | NXmx | s | | `sensor_thickness` | NXmx | m | diff --git a/rugnux/RugnuxCalibration.cpp b/rugnux/RugnuxCalibration.cpp index 39bba67d..3aa24e1f 100644 --- a/rugnux/RugnuxCalibration.cpp +++ b/rugnux/RugnuxCalibration.cpp @@ -90,14 +90,20 @@ void WritePoniFile(const std::string &path, const DiffractionExperiment &experim // pixel, i.e. it is the point of normal incidence, so it maps straight across with no correction. // GetDirectBeam_pxl() is a different quantity - where the direct beam lands - and parts from the // PONI as soon as rot1/rot2 are non-zero. + // + // The half pixel is the origin convention (see docs/DETECTOR_GEOMETRY.md): our coordinates are + // pixel-centred, so beam_x = 948 means the CENTRE of pixel 948, while pyFAI measures from the edge + // of the sensor and puts the centre of pixel i at (i + 0.5) * pixel size. Without it the pattern + // pyFAI integrates sits half a pixel off ours. + const double half_pixel_m = 0.5 * pixel_m; f << fmt::format("# Calibration done by Jungfraujoch rugnux {}\n", jfjoch_version()); f << "poni_version: 2\n"; f << "Detector: Detector\n"; f << fmt::format("Detector_config: {{\"pixel1\": {:g}, \"pixel2\": {:g}, \"max_shape\": [{}, {}]}}\n", pixel_m, pixel_m, experiment.GetYPixelsNumConv(), experiment.GetXPixelsNumConv()); f << fmt::format("Distance: {:.9g}\n", geom.GetDetectorDistance_mm() * 1e-3); - f << fmt::format("Poni1: {:.9g}\n", geom.GetBeamY_pxl() * pixel_m); - f << fmt::format("Poni2: {:.9g}\n", geom.GetBeamX_pxl() * pixel_m); + f << fmt::format("Poni1: {:.9g}\n", geom.GetBeamY_pxl() * pixel_m + half_pixel_m); + f << fmt::format("Poni2: {:.9g}\n", geom.GetBeamX_pxl() * pixel_m + half_pixel_m); // rot2 and rot3 change SIGN on the way out, and rot1 does not. pyFAI has the slow axis increasing // BOTTOM to TOP; we use the MX convention, top to bottom. The two frames therefore differ by a // reflection in y, and conjugating a rotation by a reflection gives R(n, theta) -> R(Mn, -theta). @@ -105,8 +111,8 @@ void WritePoniFile(const std::string &path, const DiffractionExperiment &experim // (about x) and rot3 (about the beam) the axis lies in the mirror plane, so only the sense reverses. // The angles mean the same thing in both frames - it is only the handedness of the frame that // differs - and the same two flips would apply on the way IN if a PONI file were ever read. - // Poni1/Poni2 need no such change: they are distances from pixel (0, 0) along each axis, which the - // direction the axis runs in does not affect. + // Poni1/Poni2 need no such change: they are distances from the corner of the sensor along each + // axis, which the direction the axis runs in does not affect. // Verified against pyFAI on a LaB6 image: written without the flip, the rings pyFAI integrates are // BROADER than with no tilt at all (peak 42 against 30, mean ring-position error 0.0045 1/A against // 0.0027); with it they sharpen to 132 and 0.0005. diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp index f287256f..4427c0db 100644 --- a/tests/CalibrationTest.cpp +++ b/tests/CalibrationTest.cpp @@ -67,8 +67,11 @@ TEST_CASE("Calibration_PoniFileAxisConvention", "[DetGeomCalib]") { const double pixel_m = geom.GetPixelSize_mm() * 1e-3; CHECK(keys["poni_version"] == "2"); - CHECK(std::stod(keys["Poni1"]) == Catch::Approx(1275.0 * pixel_m)); // slow axis = y - CHECK(std::stod(keys["Poni2"]) == Catch::Approx(1000.0 * pixel_m)); // fast axis = x + // The half pixel is the origin convention (docs/DETECTOR_GEOMETRY.md): our beam centre is + // pixel-centred, pyFAI measures from the edge of the sensor and puts the centre of pixel i at + // (i + 0.5) * pixel size. + CHECK(std::stod(keys["Poni1"]) == Catch::Approx(1275.5 * pixel_m)); // slow axis = y + CHECK(std::stod(keys["Poni2"]) == Catch::Approx(1000.5 * pixel_m)); // fast axis = x CHECK(std::stod(keys["Distance"]) == Catch::Approx(0.150)); // rot2 and rot3 are NEGATED into pyFAI's frame and rot1 is not: pyFAI's slow axis runs bottom to // top where the MX convention runs top to bottom, so the frames differ by a reflection in y. That diff --git a/tests/DiffractionGeometryTest.cpp b/tests/DiffractionGeometryTest.cpp index 306304b5..872f6963 100644 --- a/tests/DiffractionGeometryTest.cpp +++ b/tests/DiffractionGeometryTest.cpp @@ -264,17 +264,22 @@ Wavelength: 1e-10 */ // PyFAI uses nm^-1 for Q? +// The beam centre is Poni/pixel_size - 0.5 in every PONI test here: our coordinates are pixel-centred +// (0.0 is the centre of the first pixel) while pyFAI measures from the edge of the sensor and puts the +// centre of pixel i at (i + 0.5) * pixel size - see docs/DETECTOR_GEOMETRY.md. So 0.150 m / 75 um gives +// 1999.5, not 2000. With the half pixel the reference values below are reproduced to float precision; +// without it every one of them is out by 2.6e-3 nm^-1, which the old 1e-2 tolerance hid. DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV); + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV); DiffractionGeometry geom = x.GetDiffractionGeometry(); float diff_800_400 = fabs(geom.PxlToQ( 800,400)*10.0 - 6.295358803860941); float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 7.554628215027982); float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 5.73479724964891); - REQUIRE(diff_800_400 < 0.01); - REQUIRE(diff_400_800 < 0.01); - REQUIRE(diff_1300_2000 < 0.01); + REQUIRE(diff_800_400 < 1e-4); + REQUIRE(diff_400_800 < 1e-4); + REQUIRE(diff_1300_2000 < 1e-4); } TEST_CASE("DiffractionGeometry_PONI_phi","") { @@ -293,7 +298,7 @@ Wavelength: 1e-10 // PyFAI uses nm^-1 for Q? DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV); + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV); DiffractionGeometry geom = x.GetDiffractionGeometry(); float phi_2000_0 = fabs(geom.Phi_rad(2000,0) - 2 * M_PI + 1.5702959937284997); @@ -301,10 +306,10 @@ Wavelength: 1e-10 float phi_0_1000 = fabs(geom.Phi_rad(0,1000) - 3.1413425992666903); float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 2.1809518509415025); - CHECK(phi_2000_0 < 0.001); - CHECK(phi_2000_2000 < 0.001); - CHECK(phi_0_1000 < 0.001); - CHECK(phi_2000_1300 < 0.001); + CHECK(phi_2000_0 < 1e-4); + CHECK(phi_2000_2000 < 1e-4); + CHECK(phi_0_1000 < 1e-4); + CHECK(phi_2000_1300 < 1e-4); } @@ -324,7 +329,7 @@ Wavelength: 1e-10 // PyFAI uses nm^-1 for Q? DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV) + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV) .PoniRot3_rad(0.5); DiffractionGeometry geom = x.GetDiffractionGeometry(); REQUIRE(geom.GetPoniRot3_rad() == Catch::Approx(0.5f)); @@ -332,8 +337,8 @@ Wavelength: 1e-10 float phi_800_400 = fabs(geom.Phi_rad(800,400) - 3.105073518019684); float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 1.6809518509415027); - CHECK(phi_800_400 < 0.001); - CHECK(phi_2000_1300 < 0.001); + CHECK(phi_800_400 < 1e-4); + CHECK(phi_2000_1300 < 1e-4); } @@ -353,7 +358,7 @@ Wavelength: 1e-10 // PyFAI uses nm^-1 for Q? DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV) + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV) .PoniRot1_rad(0.2).PoniRot2_rad(-0.1).PoniRot3_rad(0.5); DiffractionGeometry geom = x.GetDiffractionGeometry(); @@ -364,8 +369,8 @@ Wavelength: 1e-10 float phi_800_400 = fabs(geom.Phi_rad(800,400) - 2 * M_PI + 1.4175001633470816); float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 2 * M_PI + 0.6630282166663707); - CHECK(phi_800_400 < 0.001); - CHECK(phi_2000_1300 < 0.001); + CHECK(phi_800_400 < 1e-4); + CHECK(phi_2000_1300 < 1e-4); } TEST_CASE("DiffractionGeometry_PONI_rot1","") { @@ -384,7 +389,7 @@ Wavelength: 1e-10 // PyFAI uses nm^-1 for Q? DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV); + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV); DiffractionGeometry geom = x.GetDiffractionGeometry(); geom.PoniRot1_rad(0.2); @@ -392,9 +397,9 @@ Wavelength: 1e-10 float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 5.148411999405654); float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 10.37635963741911); - CHECK(diff_800_400 < 0.01); - CHECK(diff_400_800 < 0.01); - CHECK(diff_1300_2000 < 0.01); + CHECK(diff_800_400 < 1e-4); + CHECK(diff_400_800 < 1e-4); + CHECK(diff_1300_2000 < 1e-4); } @@ -414,7 +419,7 @@ Wavelength: 1e-10 // PyFAI uses nm^-1 for Q? DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(1000).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV); + x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV); DiffractionGeometry geom = x.GetDiffractionGeometry(); geom.PoniRot1_rad(0.2).PoniRot2_rad(-0.1); @@ -422,9 +427,9 @@ Wavelength: 1e-10 float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 8.805012278158177); float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 9.363455481328781); - CHECK(diff_800_400 < 0.01); - CHECK(diff_400_800 < 0.01); - CHECK(diff_1300_2000 < 0.01); + CHECK(diff_800_400 < 1e-4); + CHECK(diff_400_800 < 1e-4); + CHECK(diff_1300_2000 < 1e-4); } TEST_CASE("DiffractionGeometry_PyFAI_Solid_angle","") { @@ -445,14 +450,14 @@ Wavelength: 1e-10 // so it is independent of the poni rotation (tilt). CalcAzIntSolidAngleCorr matches this; // the invariance is checked in DiffractionGeometry_SolidAngleCorrection_TiltInvariant. DiffractionExperiment x(DetJF4M()); - x.DetectorDistance_mm(200).BeamX_pxl(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV); + x.DetectorDistance_mm(200).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV); DiffractionGeometry geom = x.GetDiffractionGeometry(); float diff_100_100 = fabs(geom.CalcAzIntSolidAngleCorr( 100,100) - 0.4844596502755233); - CHECK(diff_100_100 < 0.0002); + CHECK(diff_100_100 < 1e-5); float diff_400_800 = fabs(geom.CalcAzIntSolidAngleCorr( 400,800)- 0.6267921080721112); - CHECK(diff_400_800 < 0.0002); + CHECK(diff_400_800 < 1e-5); } TEST_CASE("ResPhiToPxl") { diff --git a/viewer/image_viewer/JFJochDiffractionImage.cpp b/viewer/image_viewer/JFJochDiffractionImage.cpp index 7ab70443..790602b8 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.cpp +++ b/viewer/image_viewer/JFJochDiffractionImage.cpp @@ -390,7 +390,10 @@ void JFJochDiffractionImage::DrawBeamCenter() { auto geom = image->Dataset().experiment.GetDiffractionGeometry(); auto [beam_x, beam_y] = geom.GetDirectBeam_pxl(); - DrawCross(beam_x, beam_y, 25, 5, 2); + // + 0.5 as everywhere else in the overlay: our coordinates are pixel-centred, the scene's are + // pixel-cornered (pixel i covers [i, i+1)), so the cross would otherwise sit half a pixel off the + // spots and the image. + DrawCross(beam_x + 0.5f, beam_y + 0.5f, 25, 5, 2); } void JFJochDiffractionImage::DrawTopPixels() { diff --git a/viewer/windows/JFJochCalibrationResultWindow.cpp b/viewer/windows/JFJochCalibrationResultWindow.cpp index 8000c81d..98453230 100644 --- a/viewer/windows/JFJochCalibrationResultWindow.cpp +++ b/viewer/windows/JFJochCalibrationResultWindow.cpp @@ -96,9 +96,12 @@ JFJochCalibrationResultWindow::JFJochCalibrationResultWindow(const QString &titl auto *note = new QLabel(this); note->setWordWrap(true); note->setStyleSheet("color: gray;"); + // In mm these are what the .poni file carries, so they are measured pyFAI's way - from the edge of + // the sensor, which puts the centre of pixel i at (i + 0.5) * pixel size - not from the centre of + // the first pixel as the [px] rows above (docs/DETECTOR_GEOMETRY.md). note->setText(QStringLiteral("PONI x = %1 mm, PONI y = %2 mm.%3") - .arg(g.GetBeamX_pxl() * pxl_mm, 0, 'f', 4) - .arg(g.GetBeamY_pxl() * pxl_mm, 0, 'f', 4) + .arg((g.GetBeamX_pxl() + 0.5) * pxl_mm, 0, 'f', 4) + .arg((g.GetBeamY_pxl() + 0.5) * pxl_mm, 0, 'f', 4) .arg(poni_path.isEmpty() ? QString() : " Written to " + poni_path)); layout->addWidget(note); }