rugnux reads marCCD sweeps natively
A decade of deposited CCD data is archived as marCCD - what Rayonix MX-series and mar Mosaic detectors write - and rugnux could not open any of it. reader/ had two formats, NXmx/HDF5 and PILATUS miniCBF, and rugnux_cli dispatched on the one CanRead(); this adds the third. The format needs no new dependency: a marCCD file is an ordinary uncompressed TIFF whose 3072-byte instrument header sits in the gap between the TIFF header and the pixels, so libtiff - already fetched for JFJochPreview in every build mode - reads the image, and the header is a fixed-offset block of little-endian int32. Two things differ from the miniCBF path and are worth naming: * The pixel size is NOT rounded to whole micrometres. A PILATUS pixel is exactly 172 um so the existing reader can afford lround(); a MAR300 pixel is 73.242 um, and rounding it to 73 is a 0.33% scale error on every cell edge reported. * The sweep template is the last run of digits in the whole file name rather than in the stem, which covers both schemes these detectors use - a numbered stem (xtal_1_00042.mccd) and the frame number as the extension (D1.042). A CCD frame marks no untrusted pixels, so the sweep starts with nothing masked, and the format has nowhere to state the rotation axis' direction, so the run settles its sign from the data exactly as it does for a miniCBF carrying no axis table. Measured on one deposited 300-frame Rayonix MX-300 sweep, de novo with no flags: 100% indexing, the deposited point group, cell within 0.045%, 99.5% complete at multiplicity 3.4, in 26 s. The chosen sweep is confirmed against the instrument header before it is opened, so a directory of ordinary TIFFs is refused rather than read with a pixel size of zero - a unit test covers that, both naming schemes, and the geometry conversion. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_fandClaude Opus 5 committed
1 parent
cb5a2f032a
commit
a1889c45e9
9 files changed
+850
-8
No files matched your search
@@ -1,6 +1,10 @@
|
||||
# Changelog
|
||||
## 1.0.0
|
||||
|
||||
### 1.0.0-rc.171
|
||||
|
||||
* rugnux reads marCCD images natively: naming any frame of a Rayonix or mar Mosaic sweep processes the whole sweep, with no conversion step.
|
||||
|
||||
### 1.0.0-rc.170
|
||||
|
||||
* Fixed a `jfjoch_broker` crash during indexing: sorting no longer misbehaves on non-finite values, and GPU FFT indexer kernel launches are now error-checked.
|
||||
|
||||
+13
-3
@@ -1,10 +1,10 @@
|
||||
# What rugnux reads
|
||||
|
||||
Which data `rugnux` opens, before anything is typed. The short answer: an HDF5 master (NXmx or
|
||||
DECTRIS, from any facility) or a PILATUS miniCBF sweep — nothing else is read, so any other format
|
||||
has to be converted to one of these two first.
|
||||
DECTRIS, from any facility), a PILATUS miniCBF sweep or a marCCD sweep — nothing else is read, so
|
||||
any other format has to be converted to one of these three first.
|
||||
|
||||
**Input** is either an HDF5 master file or a directory of PILATUS miniCBF frames. One input is
|
||||
**Input** is an HDF5 master file, or a directory of PILATUS miniCBF or marCCD frames. One input is
|
||||
**one sweep of one crystal** — rugnux does not combine sweeps or crystals in a run; process each
|
||||
sweep to its own `_unmerged.mtz` and merge them downstream
|
||||
(see [Taking the data onward](RUGNUX_INTEGRATION.md#taking-the-data-onward)).
|
||||
@@ -26,6 +26,16 @@ sweep to its own `_unmerged.mtz` and merge them downstream
|
||||
header, including the imgCIF axis table where the header carries one (see
|
||||
[Detector geometry](DETECTOR_GEOMETRY.md)). A raw CBF carries no analysis results, so
|
||||
`--mode scale` — which re-scales the reflections stored in a `_process.h5` — does not accept one.
|
||||
* **marCCD sweep** — what Rayonix MX-series and mar Mosaic detectors write, and what a decade of
|
||||
deposited CCD data is archived as: an uncompressed TIFF with the instrument header in the gap
|
||||
before the pixels. One frame per file, read natively. Naming a frame or its directory selects the
|
||||
sweep exactly as for miniCBF, under either naming scheme these detectors use — a numbered stem
|
||||
(`xtal_1_00042.mccd`) or the frame number as the file extension (`D1.042`). The distance, beam
|
||||
centre, pixel size, wavelength and the circle that turned come from the header; the rotation
|
||||
axis' direction does not, because the format has nowhere to state it, so the run settles its sign
|
||||
from the data as it does for a miniCBF that carries no axis table. A CCD frame marks no untrusted
|
||||
pixels, so the sweep starts with nothing masked. Like a raw CBF, it carries no analysis results
|
||||
and `--mode scale` does not accept one.
|
||||
|
||||
Spots are always found by `rugnux` itself, including for the two-pass rotation first pass — the
|
||||
spot lists a dataset may already carry were found online, at the acquisition's threshold and with
|
||||
|
||||
@@ -6,6 +6,10 @@ ADD_LIBRARY(JFJochReader STATIC
|
||||
MiniCBF.h
|
||||
JFJochCBFReader.cpp
|
||||
JFJochCBFReader.h
|
||||
MarCCD.cpp
|
||||
MarCCD.h
|
||||
JFJochMarCCDReader.cpp
|
||||
JFJochMarCCDReader.h
|
||||
HDF5ImageLocator.cpp
|
||||
HDF5ImageLocator.h
|
||||
HDF5ImageSource.cpp
|
||||
@@ -22,6 +26,8 @@ ADD_LIBRARY(JFJochReader STATIC
|
||||
# NB: JFJochHttpReader (the viewer's live-broker HTTP client) lives in viewer/ and links libcurl;
|
||||
# it is deliberately NOT part of this always-built library, so the broker/writer never pull in a
|
||||
# TLS/Kerberos stack.
|
||||
# tiff: the marCCD reader's files are ordinary uncompressed TIFFs, so libtiff - which the build
|
||||
# already fetches for JFJochPreview in every mode - reads the pixels and there is no new dependency.
|
||||
TARGET_LINK_LIBRARIES(JFJochReader JFJochImageAnalysis JFJochAPI JFJochCommon JFJochZMQ JFJochLogger
|
||||
JFJochHDF5Wrappers CBORStream2FrameSerialize
|
||||
JFJochHDF5Wrappers CBORStream2FrameSerialize tiff
|
||||
${CMAKE_DL_LIBS})
|
||||
@@ -0,0 +1,190 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "JFJochMarCCDReader.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <future>
|
||||
#include <thread>
|
||||
|
||||
#include "../common/JFJochException.h"
|
||||
#include "../common/JFJochMath.h"
|
||||
#include "../common/Logger.h"
|
||||
|
||||
namespace {
|
||||
|
||||
// The base rotation axis in the internal frame (x along increasing detector column, y along
|
||||
// increasing row, z along the beam). A marCCD header names the circle that turned but never states
|
||||
// a direction for it, so this is the convention an NXmx master writes for the same instruments, and
|
||||
// a file that needs the other sign is settled from the data by the run's axis-sign rescue - the
|
||||
// same arrangement JFJochCBFReader makes for a miniCBF that states no axis table.
|
||||
const Coord ASSUMED_BASE_AXIS(-1.0f, 0.0f, 0.0f);
|
||||
|
||||
} // namespace
|
||||
|
||||
bool JFJochMarCCDReader::CanRead(const std::string &path) {
|
||||
return marccd::CanRead(path);
|
||||
}
|
||||
|
||||
void JFJochMarCCDReader::ReadFiles(const std::string &path) {
|
||||
files_ = marccd::CollectSweep(path);
|
||||
if (files_.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"No marCCD images found for " + path);
|
||||
|
||||
header0_ = marccd::ReadHeader(files_[0]);
|
||||
// A pixel size is what makes the file an IMAGE: every resolution, every scattering vector and
|
||||
// the beam centre in millimetres scale by it, and a default of 0 collapses all of them without
|
||||
// a word.
|
||||
if (!(header0_.pixel_x_m > 0.0) || !(header0_.pixel_y_m > 0.0))
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
files_[0] + " states no pixel size in its marCCD header");
|
||||
if (!(header0_.wavelength_A > 0.0))
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
files_[0] + " states no wavelength in its marCCD header");
|
||||
if (!(header0_.distance_m > 0.0))
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
files_[0] + " states no detector distance in its marCCD header");
|
||||
|
||||
dataset_ = std::make_shared<JFJochReaderDataset>();
|
||||
dataset_->experiment = default_experiment;
|
||||
|
||||
DetectorSetup detector = DetDECTRIS(header0_.nx, header0_.ny,
|
||||
header0_.detector.empty() ? "marCCD" : header0_.detector, {});
|
||||
// Not rounded to whole micrometres, as the miniCBF path can afford to be: a PILATUS pixel is
|
||||
// exactly 172 um, but these are 73.242 um, and rounding that to 73 is a 0.33% scale error on
|
||||
// every cell edge the run reports.
|
||||
detector.PixelSize_um(static_cast<float>(header0_.pixel_x_m * 1e6));
|
||||
// A CCD has no sensor thickness worth correcting for: the phosphor converts at the surface and
|
||||
// the fibre optic carries light, not X-rays, so the parallax correction a silicon sensor needs
|
||||
// does not apply. Left at zero, which is what the geometry means by "no depth".
|
||||
detector.SensorThickness_um(0);
|
||||
if (header0_.saturated_value > 0)
|
||||
detector.SaturationLimit(SaturationLimitFromValue(header0_.saturated_value));
|
||||
else
|
||||
Logger("MarCCDReader").Warning("{} states no saturated value, so no pixel will be called "
|
||||
"saturated; if this detector overloads, its strongest "
|
||||
"reflections will be integrated as if they were valid.",
|
||||
files_[0]);
|
||||
// Images are handed out as signed 32-bit whatever the file stored, so that is the depth the
|
||||
// rest of the code must see; the real overflow is the header's saturated value, set above.
|
||||
detector.BitDepthImage(32);
|
||||
detector.MinFrameTime(std::chrono::microseconds(0));
|
||||
detector.MinCountTime(std::chrono::microseconds(0));
|
||||
detector.ReadOutTime(std::chrono::nanoseconds(0));
|
||||
dataset_->experiment.Detector(detector);
|
||||
|
||||
dataset_->experiment.BeamX_pxl(static_cast<float>(header0_.beam_x_px));
|
||||
dataset_->experiment.BeamY_pxl(static_cast<float>(header0_.beam_y_px));
|
||||
dataset_->experiment.DetectorDistance_mm(static_cast<float>(header0_.distance_m * 1000.0));
|
||||
|
||||
// A detector swung out on a 2theta arm. The arm turns the detector about the sample and so
|
||||
// carries the square-on geometry with it: the header's distance stays the distance along the
|
||||
// detector normal and the beam centre stays the point of normal incidence, which is exactly
|
||||
// what the PONI convention wants, so the swing is a PONI rotation and nothing else changes.
|
||||
// The arm turns about the same axis as the spindle on the geometries these headers describe.
|
||||
if (header0_.two_theta_deg != 0.0) {
|
||||
float rot1 = 0, rot2 = 0, rot3 = 0;
|
||||
PoniAnglesFromMatrix(RotMatrix(static_cast<float>(header0_.two_theta_deg * PI / 180.0),
|
||||
ASSUMED_BASE_AXIS), rot1, rot2, rot3);
|
||||
dataset_->experiment.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3);
|
||||
}
|
||||
|
||||
dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast<float>(header0_.wavelength_A));
|
||||
dataset_->experiment.FrameTime(
|
||||
std::chrono::duration_cast<std::chrono::nanoseconds>(
|
||||
std::chrono::duration<double>(header0_.exposure_s)),
|
||||
std::chrono::duration_cast<std::chrono::nanoseconds>(
|
||||
std::chrono::duration<double>(header0_.exposure_s)));
|
||||
|
||||
// The rotation angle of every image, from its own header. Reading one costs a 4 kB read, so on
|
||||
// a sweep of several thousand frames this is worth spreading over the cores, as the CBF path
|
||||
// does for the same reason.
|
||||
std::vector<double> angles(files_.size());
|
||||
{
|
||||
const size_t nthreads = std::min<size_t>(std::max(1u, std::thread::hardware_concurrency()), 8);
|
||||
std::vector<std::future<void>> futures;
|
||||
for (size_t t = 0; t < nthreads; t++)
|
||||
futures.push_back(std::async(std::launch::async, [&, t] {
|
||||
for (size_t i = t; i < files_.size(); i += nthreads)
|
||||
angles[i] = marccd::ReadHeader(files_[i]).start_angle_deg;
|
||||
}));
|
||||
for (auto &f : futures)
|
||||
f.get();
|
||||
}
|
||||
|
||||
double increment = header0_.angle_increment_deg;
|
||||
if (files_.size() > 1) {
|
||||
// Prefer the measured step over the header's nominal one, and unwrap a sweep that passes 360.
|
||||
double d = angles[1] - angles[0];
|
||||
if (d < -180.0) d += 360.0;
|
||||
if (std::abs(d) > 1e-6) increment = d;
|
||||
}
|
||||
dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
|
||||
static_cast<float>(angles.front()),
|
||||
static_cast<float>(increment),
|
||||
ASSUMED_BASE_AXIS, {}));
|
||||
|
||||
dataset_->error_value = -1;
|
||||
dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(files_.size()));
|
||||
// A CCD frame stores no untrusted-pixel marker - every value is a real reading, and the
|
||||
// detector has no module gaps - so the sweep starts with nothing masked.
|
||||
dataset_->pixel_mask = std::make_shared<const PixelMask>(static_cast<size_t>(header0_.nx),
|
||||
static_cast<size_t>(header0_.ny));
|
||||
|
||||
SetStartMessage(dataset_);
|
||||
}
|
||||
|
||||
uint64_t JFJochMarCCDReader::GetNumberOfImages() const {
|
||||
return files_.size();
|
||||
}
|
||||
|
||||
void JFJochMarCCDReader::Close() {
|
||||
files_.clear();
|
||||
dataset_.reset();
|
||||
}
|
||||
|
||||
template <class Buffer>
|
||||
CompressedImage JFJochMarCCDReader::DecodeInto(int64_t image_number, Buffer &buffer,
|
||||
std::vector<uint8_t> &scratch) const {
|
||||
if (image_number < 0 || static_cast<size_t>(image_number) >= files_.size())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"Image number out of range");
|
||||
|
||||
const size_t npixel = static_cast<size_t>(header0_.nx) * static_cast<size_t>(header0_.ny);
|
||||
buffer.resize(npixel * sizeof(int32_t));
|
||||
|
||||
const auto h = marccd::ReadInto(files_[image_number],
|
||||
reinterpret_cast<int32_t *>(buffer.data()), npixel, scratch);
|
||||
if (h.nx != header0_.nx || h.ny != header0_.ny)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
"marCCD image size differs from the first image of the sweep");
|
||||
|
||||
return CompressedImage(buffer.data(), buffer.size(),
|
||||
static_cast<size_t>(header0_.nx), static_cast<size_t>(header0_.ny),
|
||||
CompressedImageMode::Int32, CompressionAlgorithm::NO_COMPRESSION);
|
||||
}
|
||||
|
||||
bool JFJochMarCCDReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
||||
DataMessage &message,
|
||||
std::vector<uint8_t> &buffer,
|
||||
int64_t image_number,
|
||||
bool update_dataset) {
|
||||
(void) update_dataset;
|
||||
if (!dataset)
|
||||
return false;
|
||||
|
||||
std::vector<uint8_t> scratch;
|
||||
message.image = DecodeInto(image_number, buffer, scratch);
|
||||
message.number = image_number;
|
||||
return true;
|
||||
}
|
||||
|
||||
bool JFJochMarCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
|
||||
image.image = DecodeInto(image_number, image.image_buffer, image.read_buffer);
|
||||
return true;
|
||||
}
|
||||
|
||||
std::vector<SpotToSave> JFJochMarCCDReader::ReadSpots(int64_t) const {
|
||||
return {}; // a raw marCCD file stores no analysis results
|
||||
}
|
||||
@@ -0,0 +1,57 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "JFJochReader.h"
|
||||
#include "MarCCD.h"
|
||||
|
||||
// Reads a rotation sweep straight from a directory of marCCD files - what Rayonix MX-series and mar
|
||||
// Mosaic detectors write, and what a decade of deposited CCD data is archived as - with no
|
||||
// conversion step.
|
||||
//
|
||||
// The same shape as JFJochCBFReader: the sweep's geometry comes from the first file's header, the
|
||||
// rotation angle of each image from its own, and images are decoded on demand so any number of
|
||||
// workers can read at once. A raw marCCD file carries no analysis results, so there are no spots
|
||||
// and no reflections here either.
|
||||
//
|
||||
// The one thing a CCD needs that a pixel detector does not: its pixels are unsigned 16-bit with no
|
||||
// negative marker for untrusted ones, so the mask is what the header calls saturated and nothing
|
||||
// else. A CCD's point spread and read-out noise are left to the analysis, which measures the
|
||||
// background from the image rather than assuming a detector.
|
||||
class JFJochMarCCDReader : public JFJochReader {
|
||||
std::vector<std::string> files_;
|
||||
std::shared_ptr<JFJochReaderDataset> dataset_;
|
||||
marccd::Header header0_;
|
||||
|
||||
bool LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
|
||||
DataMessage &message,
|
||||
std::vector<uint8_t> &buffer,
|
||||
int64_t image_number,
|
||||
bool update_dataset) override;
|
||||
|
||||
// Decodes one image into the caller's byte buffer and returns the image that points at it. The
|
||||
// stored 16-bit strip is read through scratch, which the caller keeps between frames.
|
||||
template <class Buffer>
|
||||
CompressedImage DecodeInto(int64_t image_number, Buffer &buffer, std::vector<uint8_t> &scratch) const;
|
||||
|
||||
public:
|
||||
~JFJochMarCCDReader() override = default;
|
||||
|
||||
// True if the path names something this reader can open: a marCCD file, or a directory holding
|
||||
// at least one. Cheap - it reads a few kB at most.
|
||||
static bool CanRead(const std::string &path);
|
||||
|
||||
// path is a directory of marCCD frames, or one frame inside the sweep to take the whole sweep
|
||||
// from.
|
||||
void ReadFiles(const std::string &path);
|
||||
|
||||
[[nodiscard]] uint64_t GetNumberOfImages() const override;
|
||||
void Close() override;
|
||||
|
||||
bool ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) override;
|
||||
[[nodiscard]] std::vector<SpotToSave> ReadSpots(int64_t image) const override;
|
||||
};
|
||||
@@ -0,0 +1,322 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "MarCCD.h"
|
||||
|
||||
#include <tiffio.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cctype>
|
||||
#include <cstring>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <map>
|
||||
#include <optional>
|
||||
#include <tuple>
|
||||
|
||||
#include "../common/JFJochException.h"
|
||||
|
||||
namespace marccd {
|
||||
|
||||
namespace {
|
||||
|
||||
// The instrument header, as marCCD writes it: 3072 bytes of little-endian int32 at fixed offsets.
|
||||
// Only the fields below are read; the rest are detector housekeeping. Offsets are from the start of
|
||||
// the block, which is where the "MMX" name sits.
|
||||
constexpr size_t HEADER_BYTES = 3072;
|
||||
constexpr size_t OFF_NAME = 4; // char[16], "MMX"
|
||||
constexpr size_t OFF_NFAST = 80;
|
||||
constexpr size_t OFF_NSLOW = 84;
|
||||
constexpr size_t OFF_DEPTH = 88; // bytes per pixel
|
||||
constexpr size_t OFF_SATURATED = 104;
|
||||
// Goniostat block: distance and beam centre, then the driven circles as start/end pairs.
|
||||
constexpr size_t OFF_XTAL_TO_DETECTOR = 640; // micrometres
|
||||
constexpr size_t OFF_BEAM_X = 644; // 1/1000 pixel, along fast
|
||||
constexpr size_t OFF_BEAM_Y = 648; // 1/1000 pixel, along slow
|
||||
constexpr size_t OFF_EXPOSURE_TIME = 656; // milliseconds
|
||||
constexpr size_t OFF_START_TWOTHETA = 668; // millidegrees, and the seven that follow it
|
||||
constexpr size_t OFF_END_TWOTHETA = 700; // the same seven, at the end of the exposure
|
||||
constexpr size_t OFF_ROTATION_RANGE = 736; // millidegrees
|
||||
constexpr size_t OFF_PIXELSIZE_X = 772; // nanometres
|
||||
constexpr size_t OFF_PIXELSIZE_Y = 776;
|
||||
constexpr size_t OFF_SOURCE_WAVELENGTH = 908; // 1e-5 angstrom
|
||||
|
||||
// The circles of the goniostat block, in the order it stores them. The scanned one is whichever
|
||||
// pair of start/end angles differs, so this table is what turns that into a name.
|
||||
constexpr const char *CIRCLE_NAMES[] = {"two_theta", "omega", "chi", "kappa",
|
||||
"phi", "delta", "gamma"};
|
||||
constexpr size_t N_CIRCLES = sizeof(CIRCLE_NAMES) / sizeof(CIRCLE_NAMES[0]);
|
||||
|
||||
int32_t I32(const std::vector<uint8_t> &h, size_t off) {
|
||||
int32_t v;
|
||||
std::memcpy(&v, h.data() + off, sizeof(v));
|
||||
return v; // the header declares little-endian byte order and no writer has ever used the other
|
||||
}
|
||||
|
||||
// The first bytes of a file: the TIFF header, the instrument header, and nothing else.
|
||||
std::vector<uint8_t> ReadPrefix(const std::string &path, size_t bytes) {
|
||||
std::ifstream f(path, std::ios::binary);
|
||||
if (!f)
|
||||
throw JFJochException(JFJochExceptionCategory::MockFileOpenError, "Cannot open " + path);
|
||||
std::vector<uint8_t> out(bytes);
|
||||
f.read(reinterpret_cast<char *>(out.data()), static_cast<std::streamsize>(bytes));
|
||||
out.resize(static_cast<size_t>(f.gcount()));
|
||||
return out;
|
||||
}
|
||||
|
||||
// Where the instrument header begins. The file states it in a private TIFF tag; every writer seen
|
||||
// puts it at 1024 and says so, but the tag is what the format defines, so it wins where present.
|
||||
uint32_t HeaderOffset(const std::vector<uint8_t> &prefix) {
|
||||
if (prefix.size() < 8)
|
||||
return DEFAULT_HEADER_OFFSET;
|
||||
const bool little = prefix[0] == 'I' && prefix[1] == 'I';
|
||||
const auto u16 = [&](size_t o) -> uint32_t {
|
||||
if (o + 2 > prefix.size()) return 0;
|
||||
return little ? (uint32_t(prefix[o]) | uint32_t(prefix[o + 1]) << 8)
|
||||
: (uint32_t(prefix[o]) << 8 | uint32_t(prefix[o + 1]));
|
||||
};
|
||||
const auto u32 = [&](size_t o) -> uint32_t {
|
||||
if (o + 4 > prefix.size()) return 0;
|
||||
return little ? (uint32_t(prefix[o]) | uint32_t(prefix[o + 1]) << 8 |
|
||||
uint32_t(prefix[o + 2]) << 16 | uint32_t(prefix[o + 3]) << 24)
|
||||
: (uint32_t(prefix[o]) << 24 | uint32_t(prefix[o + 1]) << 16 |
|
||||
uint32_t(prefix[o + 2]) << 8 | uint32_t(prefix[o + 3]));
|
||||
};
|
||||
const uint32_t ifd = u32(4);
|
||||
const uint32_t n = u16(ifd);
|
||||
for (uint32_t i = 0; i < n; i++) {
|
||||
const size_t e = ifd + 2 + static_cast<size_t>(i) * 12;
|
||||
if (e + 12 > prefix.size())
|
||||
break;
|
||||
if (u16(e) == TIFFTAG_MARCCD_HEADER_OFFSET)
|
||||
return u32(e + 8);
|
||||
}
|
||||
return DEFAULT_HEADER_OFFSET;
|
||||
}
|
||||
|
||||
// The instrument header block, or an empty vector where the file does not carry one.
|
||||
std::vector<uint8_t> InstrumentHeader(const std::string &path) {
|
||||
const uint32_t off = HeaderOffset(ReadPrefix(path, DEFAULT_HEADER_OFFSET));
|
||||
auto prefix = ReadPrefix(path, off + HEADER_BYTES);
|
||||
if (prefix.size() < off + HEADER_BYTES)
|
||||
return {};
|
||||
std::vector<uint8_t> h(prefix.begin() + off, prefix.begin() + off + HEADER_BYTES);
|
||||
if (std::memcmp(h.data() + OFF_NAME, "MMX", 4) != 0)
|
||||
return {};
|
||||
return h;
|
||||
}
|
||||
|
||||
Header Parse(const std::vector<uint8_t> &h, const std::string &path) {
|
||||
Header out;
|
||||
out.nx = I32(h, OFF_NFAST);
|
||||
out.ny = I32(h, OFF_NSLOW);
|
||||
if (out.nx <= 0 || out.ny <= 0)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
path + ": marCCD header states no image size");
|
||||
if (I32(h, OFF_DEPTH) != 2)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
path + ": only 16-bit marCCD images are supported");
|
||||
|
||||
out.pixel_x_m = I32(h, OFF_PIXELSIZE_X) * 1e-9;
|
||||
out.pixel_y_m = I32(h, OFF_PIXELSIZE_Y) * 1e-9;
|
||||
out.distance_m = I32(h, OFF_XTAL_TO_DETECTOR) * 1e-6;
|
||||
// The beam centre is stored in thousandths of a pixel, x along the fast direction and y along
|
||||
// the slow one, which is the same order the internal frame counts columns and rows in.
|
||||
out.beam_x_px = I32(h, OFF_BEAM_X) * 1e-3;
|
||||
out.beam_y_px = I32(h, OFF_BEAM_Y) * 1e-3;
|
||||
out.wavelength_A = I32(h, OFF_SOURCE_WAVELENGTH) * 1e-5;
|
||||
out.two_theta_deg = I32(h, OFF_START_TWOTHETA) * 1e-3;
|
||||
out.exposure_s = I32(h, OFF_EXPOSURE_TIME) * 1e-3;
|
||||
out.saturated_value = I32(h, OFF_SATURATED);
|
||||
|
||||
// Which circle moved. The header also carries a rotation_axis index, but the start/end pair
|
||||
// that differs is the same answer read off the angles themselves, and it stays right on a file
|
||||
// whose index field was never filled in.
|
||||
out.angle_increment_deg = I32(h, OFF_ROTATION_RANGE) * 1e-3;
|
||||
for (size_t i = 0; i < N_CIRCLES; i++) {
|
||||
const double start = I32(h, OFF_START_TWOTHETA + 4 * i) * 1e-3;
|
||||
const double end = I32(h, OFF_END_TWOTHETA + 4 * i) * 1e-3;
|
||||
if (start != end) {
|
||||
out.axis_name = CIRCLE_NAMES[i];
|
||||
out.start_angle_deg = start;
|
||||
if (out.angle_increment_deg == 0.0)
|
||||
out.angle_increment_deg = end - start;
|
||||
break;
|
||||
}
|
||||
}
|
||||
return out;
|
||||
}
|
||||
|
||||
// The extensions a marCCD frame is written under. A frame whose extension is the frame number
|
||||
// (D1.001) is named by the mar software itself and is just as common as a named one.
|
||||
bool PlausibleExtension(const std::filesystem::path &p) {
|
||||
std::string ext = p.extension().string();
|
||||
if (ext.empty())
|
||||
return false;
|
||||
ext.erase(ext.begin()); // drop the dot
|
||||
if (std::all_of(ext.begin(), ext.end(), [](unsigned char c) { return std::isdigit(c); }))
|
||||
return true;
|
||||
std::transform(ext.begin(), ext.end(), ext.begin(),
|
||||
[](unsigned char c) { return std::tolower(c); });
|
||||
return ext == "mccd" || ext == "img" || ext == "tif" || ext == "tiff" || ext == "marccd";
|
||||
}
|
||||
|
||||
// The sweep a file belongs to, as a template: everything before the LAST run of digits in the whole
|
||||
// file name, the length of that run, and whatever follows it. That one rule covers both naming
|
||||
// schemes in use - "xtal_1_00042.mccd" -> {"xtal_1_", 5, ".mccd"} and "D1.042" -> {"D1.", 3, ""} -
|
||||
// because in both the frame number is the last number in the name.
|
||||
struct Template {
|
||||
std::string prefix;
|
||||
size_t digits = 0;
|
||||
std::string suffix;
|
||||
|
||||
bool operator<(const Template &o) const {
|
||||
return std::tie(prefix, digits, suffix) < std::tie(o.prefix, o.digits, o.suffix);
|
||||
}
|
||||
bool Matches(const std::string &name) const {
|
||||
if (name.size() != prefix.size() + digits + suffix.size())
|
||||
return false;
|
||||
if (name.compare(0, prefix.size(), prefix) != 0)
|
||||
return false;
|
||||
if (name.compare(prefix.size() + digits, suffix.size(), suffix) != 0)
|
||||
return false;
|
||||
for (size_t i = 0; i < digits; i++)
|
||||
if (!std::isdigit(static_cast<unsigned char>(name[prefix.size() + i])))
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
std::optional<Template> TemplateOf(const std::string &name) {
|
||||
size_t end = name.size();
|
||||
while (end > 0 && !std::isdigit(static_cast<unsigned char>(name[end - 1])))
|
||||
end--;
|
||||
if (end == 0)
|
||||
return {}; // no number anywhere: not part of a sweep
|
||||
size_t start = end;
|
||||
while (start > 0 && std::isdigit(static_cast<unsigned char>(name[start - 1])))
|
||||
start--;
|
||||
return Template{name.substr(0, start), end - start, name.substr(end)};
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
std::vector<std::string> CollectSweep(const std::string &path) {
|
||||
std::filesystem::path p(path);
|
||||
std::error_code ec;
|
||||
const bool is_dir = std::filesystem::is_directory(p, ec);
|
||||
std::filesystem::path dir = is_dir ? p : p.parent_path();
|
||||
if (dir.empty())
|
||||
dir = ".";
|
||||
|
||||
// Naming a frame selects ITS sweep; naming a directory selects the sweep with the most frames
|
||||
// in it, which is the one a user pointing at a data directory means.
|
||||
std::optional<Template> want;
|
||||
if (!is_dir)
|
||||
want = TemplateOf(p.filename().string());
|
||||
|
||||
std::map<Template, std::vector<std::string>> sweeps;
|
||||
for (const auto &e : std::filesystem::directory_iterator(dir, ec)) {
|
||||
if (!e.is_regular_file() || !PlausibleExtension(e.path()))
|
||||
continue;
|
||||
const std::string name = e.path().filename().string();
|
||||
const auto t = TemplateOf(name);
|
||||
if (!t.has_value())
|
||||
continue;
|
||||
if (want.has_value() && !want->Matches(name))
|
||||
continue;
|
||||
sweeps[*t].push_back(e.path().string());
|
||||
}
|
||||
|
||||
std::vector<std::string> out;
|
||||
for (auto &[key, files] : sweeps)
|
||||
if (files.size() > out.size())
|
||||
out = std::move(files);
|
||||
|
||||
// The frame number is zero-padded in every mar naming scheme seen, so within one template a
|
||||
// plain sort is the collection order.
|
||||
std::sort(out.begin(), out.end());
|
||||
// A directory of TIFFs from anything else matches the name rule just as well, so the chosen
|
||||
// sweep is confirmed against the one thing only a marCCD file has.
|
||||
if (!out.empty() && InstrumentHeader(out.front()).empty())
|
||||
return {};
|
||||
return out;
|
||||
}
|
||||
|
||||
bool CanRead(const std::string &path) {
|
||||
std::error_code ec;
|
||||
if (std::filesystem::is_directory(path, ec))
|
||||
return !CollectSweep(path).empty();
|
||||
if (!PlausibleExtension(std::filesystem::path(path)))
|
||||
return false;
|
||||
try {
|
||||
return !InstrumentHeader(path).empty();
|
||||
} catch (const JFJochException &) {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
Header ReadHeader(const std::string &path) {
|
||||
const auto h = InstrumentHeader(path);
|
||||
if (h.empty())
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
path + " carries no marCCD instrument header");
|
||||
return Parse(h, path);
|
||||
}
|
||||
|
||||
Header ReadInto(const std::string &path, int32_t *out, size_t capacity, std::vector<uint8_t> &scratch) {
|
||||
const Header header = ReadHeader(path);
|
||||
const size_t npixel = static_cast<size_t>(header.nx) * static_cast<size_t>(header.ny);
|
||||
if (npixel > capacity)
|
||||
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds,
|
||||
path + ": image is larger than the buffer given for it");
|
||||
|
||||
TIFF *tiff = TIFFOpen(path.c_str(), "r");
|
||||
if (tiff == nullptr)
|
||||
throw JFJochException(JFJochExceptionCategory::TIFFGeneratorError, "Cannot open TIFF " + path);
|
||||
struct Closer {
|
||||
TIFF *t;
|
||||
~Closer() { TIFFClose(t); }
|
||||
} closer{tiff};
|
||||
|
||||
uint32_t width = 0, length = 0;
|
||||
uint16_t bits = 0, samples = 1;
|
||||
TIFFGetField(tiff, TIFFTAG_IMAGEWIDTH, &width);
|
||||
TIFFGetField(tiff, TIFFTAG_IMAGELENGTH, &length);
|
||||
TIFFGetField(tiff, TIFFTAG_BITSPERSAMPLE, &bits);
|
||||
TIFFGetFieldDefaulted(tiff, TIFFTAG_SAMPLESPERPIXEL, &samples);
|
||||
if (width != header.nx || length != header.ny)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
path + ": TIFF image size differs from the marCCD header");
|
||||
if (bits != 16 || samples != 1)
|
||||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||||
path + ": only single-channel 16-bit marCCD images are supported");
|
||||
|
||||
const size_t scanline_bytes = static_cast<size_t>(width) * sizeof(uint16_t);
|
||||
if (scanline_bytes != static_cast<size_t>(TIFFScanlineSize(tiff)))
|
||||
throw JFJochException(JFJochExceptionCategory::TIFFGeneratorError,
|
||||
path + ": TIFFScanlineSize mismatch");
|
||||
scratch.resize(scanline_bytes * length);
|
||||
for (uint32_t row = 0; row < length; row++) {
|
||||
if (TIFFReadScanline(tiff, scratch.data() + static_cast<size_t>(row) * scanline_bytes,
|
||||
row, 0) < 0)
|
||||
throw JFJochException(JFJochExceptionCategory::TIFFGeneratorError,
|
||||
path + ": TIFFReadScanline error");
|
||||
}
|
||||
|
||||
// The stored pixels are unsigned 16-bit; the rest of the code works in signed 32-bit, so they
|
||||
// are widened here and nothing downstream has to know what the file held.
|
||||
const auto *src = reinterpret_cast<const uint16_t *>(scratch.data());
|
||||
for (size_t i = 0; i < npixel; i++)
|
||||
out[i] = static_cast<int32_t>(src[i]);
|
||||
return header;
|
||||
}
|
||||
|
||||
Header Read(const std::string &path, std::vector<int32_t> &out) {
|
||||
const Header header = ReadHeader(path);
|
||||
out.resize(static_cast<size_t>(header.nx) * static_cast<size_t>(header.ny));
|
||||
std::vector<uint8_t> scratch;
|
||||
return ReadInto(path, out.data(), out.size(), scratch);
|
||||
}
|
||||
|
||||
} // namespace marccd
|
||||
@@ -0,0 +1,66 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <cstdint>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
// marCCD: an ordinary uncompressed TIFF whose 3072-byte instrument header sits in the gap between
|
||||
// the TIFF header and the pixels. The image is read through libtiff, which the build already
|
||||
// carries; the instrument header is a fixed-offset block of little-endian int32, so it is read
|
||||
// directly. Rayonix MX-series and mar Mosaic detectors both write it, under half a dozen file
|
||||
// extensions (.mccd, .img, .tif, or a bare frame number).
|
||||
//
|
||||
// Every value the header stores is an integer times a fixed scale: lengths in micrometres or
|
||||
// nanometres, angles in millidegrees, the wavelength in 1e-5 A. They are converted here, so the
|
||||
// rest of the code sees metres, degrees and angstroms like every other reader.
|
||||
namespace marccd {
|
||||
|
||||
// Byte offset of the instrument header when the file does not say. Every writer seen puts it here
|
||||
// and states it in TIFF tag 34710 as well.
|
||||
inline constexpr uint32_t DEFAULT_HEADER_OFFSET = 1024;
|
||||
|
||||
// The private TIFF tag in which marCCD states where its own header begins.
|
||||
inline constexpr uint16_t TIFFTAG_MARCCD_HEADER_OFFSET = 34710;
|
||||
|
||||
struct Header {
|
||||
std::string detector; // from the TIFF Model/Software tags where they say anything
|
||||
int64_t nx = 0; // fast dimension (columns)
|
||||
int64_t ny = 0; // slow dimension (rows)
|
||||
double pixel_x_m = 0;
|
||||
double pixel_y_m = 0;
|
||||
double distance_m = 0;
|
||||
double beam_x_px = 0;
|
||||
double beam_y_px = 0;
|
||||
double wavelength_A = 0;
|
||||
double start_angle_deg = 0;
|
||||
double angle_increment_deg = 0;
|
||||
double two_theta_deg = 0;
|
||||
double exposure_s = 0;
|
||||
int64_t saturated_value = 0;
|
||||
std::string axis_name = "phi"; // the circle that moved between start and end
|
||||
};
|
||||
|
||||
// True if the path names something this reader can open: a marCCD file, or a directory holding at
|
||||
// least one. Reads a few kB at most.
|
||||
bool CanRead(const std::string &path);
|
||||
|
||||
// The files of one sweep, in collection order. path is a directory, or one frame inside the sweep
|
||||
// to take the whole sweep from. Empty where there is none.
|
||||
std::vector<std::string> CollectSweep(const std::string &path);
|
||||
|
||||
// Header only - reads the first few kB of the file. Cheap enough to call per frame for the
|
||||
// rotation angles.
|
||||
Header ReadHeader(const std::string &path);
|
||||
|
||||
// Header + pixels of one file. out is resized to nx * ny.
|
||||
Header Read(const std::string &path, std::vector<int32_t> &out);
|
||||
|
||||
// The same, widening the stored uint16 into memory the caller already has. scratch carries the
|
||||
// 16-bit strip between frames, so a worker walking a sweep allocates once. Throws if the image
|
||||
// does not fit in capacity pixels.
|
||||
Header ReadInto(const std::string &path, int32_t *out, size_t capacity, std::vector<uint8_t> &scratch);
|
||||
|
||||
} // namespace marccd
|
||||
+13
-4
@@ -19,6 +19,7 @@
|
||||
|
||||
#include "../reader/JFJochHDF5Reader.h"
|
||||
#include "../reader/JFJochCBFReader.h"
|
||||
#include "../reader/JFJochMarCCDReader.h"
|
||||
#include "../common/Logger.h"
|
||||
#include "../common/Definitions.h"
|
||||
#include "../common/DiffractionExperiment.h"
|
||||
@@ -1379,18 +1380,26 @@ static int RunRugnux(int argc, char **argv) {
|
||||
// and the beam-stop pre-scan instead of landing inside the first pass.
|
||||
auto cufft_warmup = std::async(std::launch::async, WarmUpCuFFT);
|
||||
|
||||
// 1. Read Input File. Either an HDF5 master, or one PILATUS miniCBF frame - naming any frame of a
|
||||
// sweep reads the whole sweep, which is how these are archived and how a user thinks of them.
|
||||
// 1. Read Input File. An HDF5 master, one PILATUS miniCBF frame, or one marCCD frame - naming
|
||||
// any frame of a sweep reads the whole sweep, which is how these are archived and how a user
|
||||
// thinks of them.
|
||||
JFJochHDF5Reader hdf5_reader;
|
||||
JFJochCBFReader cbf_reader;
|
||||
JFJochMarCCDReader marccd_reader;
|
||||
JFJochReader *reader_ptr = nullptr;
|
||||
const bool input_is_cbf = JFJochCBFReader::CanRead(input_file);
|
||||
const bool input_is_marccd = !input_is_cbf && JFJochMarCCDReader::CanRead(input_file);
|
||||
try {
|
||||
if (input_is_cbf) {
|
||||
cbf_reader.ReadFiles(input_file);
|
||||
reader_ptr = &cbf_reader;
|
||||
logger.Info("Read {} CBF images of the sweep containing {}",
|
||||
cbf_reader.GetNumberOfImages(), input_file);
|
||||
} else if (input_is_marccd) {
|
||||
marccd_reader.ReadFiles(input_file);
|
||||
reader_ptr = &marccd_reader;
|
||||
logger.Info("Read {} marCCD images of the sweep containing {}",
|
||||
marccd_reader.GetNumberOfImages(), input_file);
|
||||
} else {
|
||||
hdf5_reader.ReadFile(input_file);
|
||||
reader_ptr = &hdf5_reader;
|
||||
@@ -1483,9 +1492,9 @@ static int RunRugnux(int argc, char **argv) {
|
||||
"runs no space-group search. Use --mode mx, which builds the P1 cross-check "
|
||||
"merge the ledger is folded from.");
|
||||
// Re-scaling reads reflections a previous run integrated, which only a _process.h5 holds.
|
||||
if (input_is_cbf) {
|
||||
if (input_is_cbf || input_is_marccd) {
|
||||
logger.Error("--mode scale needs the integrated reflections in a _process.h5; "
|
||||
"a CBF sweep holds none. Run --mode mx on it first.");
|
||||
"a raw image sweep holds none. Run --mode mx on it first.");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
const auto total_images = static_cast<int>(reader.GetNumberOfImages());
|
||||
|
||||
@@ -8,9 +8,11 @@
|
||||
#include "../writer/FileWriter.h"
|
||||
#include "../reader/JFJochHDF5Reader.h"
|
||||
#include "../reader/JFJochCBFReader.h"
|
||||
#include "../reader/JFJochMarCCDReader.h"
|
||||
#include "../reader/MiniCBF.h"
|
||||
#include "../compression/JFJochCompressor.h"
|
||||
|
||||
#include <cstring>
|
||||
#include <fstream>
|
||||
#include <future>
|
||||
#include <iomanip>
|
||||
@@ -4113,3 +4115,179 @@ TEST_CASE("JFJochCBFReader_incomplete_header_is_refused_not_misread", "[HDF5][Fu
|
||||
RemoveMiniCBFSweep("cbfnocut");
|
||||
}
|
||||
}
|
||||
|
||||
// marCCD: a TIFF whose instrument header sits in the gap between the TIFF header and the pixels.
|
||||
// The fixtures below are written byte for byte rather than through libtiff, because the layout IS
|
||||
// what the reader has to cope with - the private tag that states where the instrument header
|
||||
// begins, and pixels that start well past the end of the IFD.
|
||||
namespace {
|
||||
constexpr size_t MARCCD_HEADER_OFF = 1024;
|
||||
constexpr size_t MARCCD_DATA_OFF = 4096;
|
||||
|
||||
void PutU16(std::vector<uint8_t> &b, size_t off, uint16_t v) {
|
||||
b[off] = v & 0xff; b[off + 1] = (v >> 8) & 0xff;
|
||||
}
|
||||
void PutU32(std::vector<uint8_t> &b, size_t off, uint32_t v) {
|
||||
for (int i = 0; i < 4; i++) b[off + i] = (v >> (8 * i)) & 0xff;
|
||||
}
|
||||
void PutI32(std::vector<uint8_t> &b, size_t off, int32_t v) {
|
||||
PutU32(b, off, static_cast<uint32_t>(v));
|
||||
}
|
||||
|
||||
struct MarCCDFields {
|
||||
int32_t start_phi_mdeg = 0;
|
||||
int32_t end_phi_mdeg = 200;
|
||||
int32_t rotation_range_mdeg = 200;
|
||||
int32_t two_theta_mdeg = 0;
|
||||
int32_t distance_um = 300000;
|
||||
int32_t beam_x_mpx = 12000; // thousandths of a pixel
|
||||
int32_t beam_y_mpx = 8000;
|
||||
int32_t pixel_nm = 73242;
|
||||
int32_t wavelength_e5A = 97872;
|
||||
int32_t saturated = 65535;
|
||||
int32_t depth = 2;
|
||||
bool write_instrument_header = true;
|
||||
};
|
||||
|
||||
// One frame. The image is all zeros: the reader is being asked about geometry, and a pixel
|
||||
// pattern would say nothing extra about that.
|
||||
void WriteMarCCDFrame(const std::string &name, int64_t nx, int64_t ny, const MarCCDFields &f) {
|
||||
const size_t npixel = static_cast<size_t>(nx) * static_cast<size_t>(ny);
|
||||
std::vector<uint8_t> file(MARCCD_DATA_OFF + npixel * 2, 0);
|
||||
|
||||
// TIFF header: little-endian, IFD at byte 8.
|
||||
file[0] = 'I'; file[1] = 'I'; PutU16(file, 2, 42); PutU32(file, 4, 8);
|
||||
struct Entry { uint16_t tag; uint16_t type; uint32_t count; uint32_t value; };
|
||||
const std::vector<Entry> ifd = {
|
||||
{256, 4, 1, static_cast<uint32_t>(nx)}, // ImageWidth
|
||||
{257, 4, 1, static_cast<uint32_t>(ny)}, // ImageLength
|
||||
{258, 3, 1, 16}, // BitsPerSample
|
||||
{259, 3, 1, 1}, // Compression: none
|
||||
{262, 3, 1, 1}, // Photometric: min-is-black
|
||||
{273, 4, 1, MARCCD_DATA_OFF}, // StripOffsets
|
||||
{277, 3, 1, 1}, // SamplesPerPixel
|
||||
{278, 4, 1, static_cast<uint32_t>(ny)}, // RowsPerStrip: one strip
|
||||
{279, 4, 1, static_cast<uint32_t>(npixel * 2)}, // StripByteCounts
|
||||
{34710, 4, 1, MARCCD_HEADER_OFF}, // where the instrument header begins
|
||||
};
|
||||
PutU16(file, 8, static_cast<uint16_t>(ifd.size()));
|
||||
size_t e = 10;
|
||||
for (const auto &t : ifd) {
|
||||
PutU16(file, e, t.tag); PutU16(file, e + 2, t.type); PutU32(file, e + 4, t.count);
|
||||
// A SHORT value sits in the first two bytes of the value field, a LONG fills it.
|
||||
if (t.type == 3) PutU16(file, e + 8, static_cast<uint16_t>(t.value));
|
||||
else PutU32(file, e + 8, t.value);
|
||||
e += 12;
|
||||
}
|
||||
PutU32(file, e, 0); // no next IFD
|
||||
|
||||
if (f.write_instrument_header) {
|
||||
const size_t h = MARCCD_HEADER_OFF;
|
||||
std::memcpy(file.data() + h + 4, "MMX", 3);
|
||||
PutI32(file, h + 80, static_cast<int32_t>(nx)); // nfast
|
||||
PutI32(file, h + 84, static_cast<int32_t>(ny)); // nslow
|
||||
PutI32(file, h + 88, f.depth);
|
||||
PutI32(file, h + 104, f.saturated);
|
||||
PutI32(file, h + 640, f.distance_um);
|
||||
PutI32(file, h + 644, f.beam_x_mpx);
|
||||
PutI32(file, h + 648, f.beam_y_mpx);
|
||||
PutI32(file, h + 656, 1000); // exposure, ms
|
||||
PutI32(file, h + 668, f.two_theta_mdeg); // start_twotheta
|
||||
PutI32(file, h + 684, f.start_phi_mdeg); // start_phi
|
||||
PutI32(file, h + 700, f.two_theta_mdeg); // end_twotheta
|
||||
PutI32(file, h + 716, f.end_phi_mdeg); // end_phi
|
||||
PutI32(file, h + 736, f.rotation_range_mdeg);
|
||||
PutI32(file, h + 772, f.pixel_nm);
|
||||
PutI32(file, h + 776, f.pixel_nm);
|
||||
PutI32(file, h + 908, f.wavelength_e5A);
|
||||
}
|
||||
|
||||
std::ofstream out(name, std::ios::binary);
|
||||
out.write(reinterpret_cast<const char *>(file.data()),
|
||||
static_cast<std::streamsize>(file.size()));
|
||||
}
|
||||
|
||||
// Two frames 0.2 deg apart, under whatever naming scheme the caller asks for.
|
||||
void WriteMarCCDSweep(const std::string &prefix, const std::string &suffix,
|
||||
int64_t nx, int64_t ny, MarCCDFields f = {}) {
|
||||
for (int frame = 1; frame <= 2; frame++) {
|
||||
f.start_phi_mdeg = (frame - 1) * 200;
|
||||
f.end_phi_mdeg = frame * 200;
|
||||
std::ostringstream name;
|
||||
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
|
||||
WriteMarCCDFrame(name.str(), nx, ny, f);
|
||||
}
|
||||
}
|
||||
|
||||
void RemoveMarCCDSweep(const std::string &prefix, const std::string &suffix) {
|
||||
for (int frame = 1; frame <= 2; frame++) {
|
||||
std::ostringstream name;
|
||||
name << prefix << std::setfill('0') << std::setw(3) << frame << suffix;
|
||||
remove(name.str().c_str());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("JFJochMarCCDReader_Geometry", "[HDF5][Full]") {
|
||||
const int64_t nx = 24, ny = 16;
|
||||
|
||||
SECTION("the header's geometry reaches the experiment in the units the rest of the code uses") {
|
||||
WriteMarCCDSweep("marccd_", ".mccd", nx, ny);
|
||||
REQUIRE(JFJochMarCCDReader::CanRead("marccd_001.mccd"));
|
||||
JFJochMarCCDReader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFiles("marccd_001.mccd"));
|
||||
CHECK(reader.GetNumberOfImages() == 2);
|
||||
|
||||
const auto x = reader.GetDataset()->experiment;
|
||||
CHECK(x.GetDetectorDistance_mm() == Catch::Approx(300.0));
|
||||
CHECK(x.GetBeamX_pxl() == Catch::Approx(12.0));
|
||||
CHECK(x.GetBeamY_pxl() == Catch::Approx(8.0));
|
||||
CHECK(x.GetPixelSize_mm() == Catch::Approx(0.073242).epsilon(1e-4));
|
||||
CHECK(x.GetWavelength_A() == Catch::Approx(0.97872).epsilon(1e-5));
|
||||
REQUIRE(x.GetGoniometer().has_value());
|
||||
// The circle that moved between start and end is the scanned one, and the step comes from
|
||||
// the two frames rather than from the header's nominal range.
|
||||
CHECK(x.GetGoniometer()->GetName() == "phi");
|
||||
CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2));
|
||||
CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(0.0));
|
||||
reader.Close();
|
||||
RemoveMarCCDSweep("marccd_", ".mccd");
|
||||
}
|
||||
|
||||
SECTION("a frame number as the file extension names a sweep just as well") {
|
||||
// What the mar software writes by default: D1.001, D1.002, ... The frame number is the last
|
||||
// number in the name either way, which is the one rule the sweep template uses.
|
||||
WriteMarCCDSweep("marccdnum.", "", nx, ny);
|
||||
REQUIRE(JFJochMarCCDReader::CanRead("marccdnum.001"));
|
||||
JFJochMarCCDReader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFiles("marccdnum.001"));
|
||||
CHECK(reader.GetNumberOfImages() == 2);
|
||||
reader.Close();
|
||||
RemoveMarCCDSweep("marccdnum.", "");
|
||||
}
|
||||
|
||||
SECTION("the pixels come back as the int32 the rest of the code reads") {
|
||||
WriteMarCCDSweep("marccdpix_", ".mccd", nx, ny);
|
||||
JFJochMarCCDReader reader;
|
||||
REQUIRE_NOTHROW(reader.ReadFiles("marccdpix_001.mccd"));
|
||||
JFJochReaderRawImage image;
|
||||
REQUIRE(reader.ReadRawImage(0, image));
|
||||
CHECK(image.image.GetWidth() == nx);
|
||||
CHECK(image.image.GetHeight() == ny);
|
||||
CHECK(image.image.GetMode() == CompressedImageMode::Int32);
|
||||
reader.Close();
|
||||
RemoveMarCCDSweep("marccdpix_", ".mccd");
|
||||
}
|
||||
|
||||
SECTION("a TIFF that is not marCCD is refused, not misread") {
|
||||
// The name rule alone matches any numbered TIFF in the directory, so a file with no
|
||||
// instrument header must not be opened with a pixel size and a distance of zero.
|
||||
MarCCDFields f;
|
||||
f.write_instrument_header = false;
|
||||
WriteMarCCDFrame("marccdplain_001.tif", nx, ny, f);
|
||||
CHECK_FALSE(JFJochMarCCDReader::CanRead("marccdplain_001.tif"));
|
||||
JFJochMarCCDReader reader;
|
||||
CHECK_THROWS(reader.ReadFiles("marccdplain_001.tif"));
|
||||
remove("marccdplain_001.tif");
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user