Jungfraujoch/reader/JFJochMarCCDReader.cpp

// 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"
#include "SweepLayout.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);
    // The marCCD field is the overload MARKER itself - the value a saturated pixel is written as -
    // not the last valid count an NXmx saturation_value states, so it is the exclusive limit
    // already and must not be raised by one. It is typically 65535 in a 16-bit image: through
    // SaturationLimitFromValue that becomes 65536, which no pixel can reach, and nothing is ever
    // called saturated. A detector block stuck at the marker then integrates as if it were signal.
    if (header0_.saturated_value > 0)
        detector.SaturationLimit(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));
    // Only when the header states a sane one. The exposure field is not always filled in: one
    // deposited sweep carries -2093438692 there, which is not a time at all, and passing it on
    // refuses the whole dataset over a number that affects no geometry and no result. A day is a
    // generous upper bound for a single frame.
    if (header0_.exposure_s > 0.0 && header0_.exposure_s < 86400.0)
        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)));
    else
        Logger("MarCCDReader").Warning("{} states an implausible exposure time ({} s); the "
                                     "default frame time is kept. Nothing in the geometry or the "
                                     "merge depends on it.", files_[0], header0_.exposure_s);

    // Where every image sits on the spindle, and how the instrument stood, 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<sweep::Frame> frames(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) {
                    const auto h = marccd::ReadHeader(files_[i]);
                    frames[i] = {files_[i], h.start_angle_deg, h.angle_increment_deg, h.distance_m,
                                 h.beam_x_px, h.beam_y_px, h.wavelength_A};
                }
            }));
        for (auto &f : futures)
            f.get();
    }

    // The sweep the headers describe, which is not always the files laid out end to end: a deposited
    // series can be missing frames, and those are gaps in the rotation rather than images to close
    // up. A folder of screening shots taken at scattered angles is refused here by name instead of
    // failing later as a lattice nobody can explain.
    const auto layout = sweep::Place(frames, "MarCCDReader");
    files_ = layout.files;
    dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
                                                   static_cast<float>(layout.start_deg),
                                                   static_cast<float>(layout.increment_deg),
                                                   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");
    if (files_[image_number].empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No image at this point of the sweep");

    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;

    if (!HasImage(image_number))
        return false;

    std::vector<uint8_t> scratch;
    message.image = DecodeInto(image_number, buffer, scratch);
    message.number = image_number;
    return true;
}

// A slot the series has no file for is a missing image, not an error: every image loop in the
// pipeline already treats "nothing to read" as a frame to pass over, which is exactly what a gap in
// a deposited sweep is.
bool JFJochMarCCDReader::HasImage(int64_t image_number) const {
    return image_number >= 0 && static_cast<size_t>(image_number) < files_.size()
           && !files_[image_number].empty();
}

bool JFJochMarCCDReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
    if (!HasImage(image_number))
        return false;
    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
}