Jungfraujoch/reader/JFJochSMVReader.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "JFJochSMVReader.h"

#include <cmath>

#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 SMV 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);

// A d*TREK vector in the internal frame. d*TREK states its vectors in the imgCIF laboratory frame
// (z from the sample to the source, y up), which differs from the internal one by a half turn about
// x - the same relation JFJochCBFReader uses for a miniCBF axis table.
Coord ImgCIFToInternal(const std::array<double, 3> &v) {
    return {static_cast<float>(v[0]), static_cast<float>(-v[1]), static_cast<float>(-v[2])};
}

}  // namespace

bool JFJochSMVReader::CanRead(const std::string &path) {
    return smv::CanRead(path);
}

void JFJochSMVReader::ReadFiles(const std::string &path) {
    files_ = smv::CollectSweep(path);
    if (files_.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No SMV images found for " + path);

    header0_ = smv::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 SMV header");
    if (!(header0_.wavelength_A > 0.0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files_[0] + " states no wavelength in its SMV header");
    if (!(header0_.distance_m > 0.0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files_[0] + " states no detector distance in its SMV header");

    dataset_ = std::make_shared<JFJochReaderDataset>();
    dataset_->experiment = default_experiment;

    DetectorSetup detector = DetDECTRIS(header0_.nx, header0_.ny,
                                        header0_.detector.empty() ? "SMV" : 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);
    // CCD_IMAGE_SATURATION is the value a saturated pixel carries, as the marCCD field is, so it is
    // the exclusive limit itself. Where the header omits it, the container's own maximum is all
    // there is to go on: it has to be said here rather than left to the fallback, because the
    // pixels are handed out as 32-bit below and the fallback would then be INT32_MAX, which no
    // 16-bit image can reach - nothing would ever be called saturated.
    const int64_t container = (1LL << (8 * header0_.bytes_per_pixel)) - 1;
    detector.SaturationLimit(header0_.saturation > 0 ? header0_.saturation : container);
    if (header0_.saturation <= 0)
        Logger("SMVReader").Warning("{}: the header states no CCD_IMAGE_SATURATION, so saturation is "
                                    "judged on the container alone; a detector that overloads below "
                                    "{} will have its strongest reflections integrated as if they "
                                    "were valid.", files_[0], container);
    // Images are handed out as signed 32-bit whatever the file stored, so that is the depth the
    // rest of the code must see.
    detector.BitDepthImage(32);
    // How the stored image sits in the detector plane, where the header states it. A Saturn's image
    // is mirrored and turned a quarter relative to the internal convention (fast x slow points at the
    // source, and the spindle runs along the rows) - the mirror is what a rotation cannot express,
    // and without it the hand of every structure would be inverted.
    if (header0_.fast_direction.has_value() && header0_.slow_direction.has_value()) {
        const auto orientation = DetectorOrientation::Match(ImgCIFToInternal(*header0_.fast_direction),
                                                            ImgCIFToInternal(*header0_.slow_direction));
        if (!orientation.has_value())
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  files_[0] + ": the image directions its header states are not a "
                                              "square-on layout this reader supports");
        detector.ImageOrientation(orientation.value());
    }
    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);
    }
    // d*TREK states every detector circle with its vector instead. They compose outermost first, as
    // they carry each other: the first named circle carries all the ones after it.
    if (!header0_.detector_circles.empty()) {
        RotMatrix r;
        for (const auto &c: header0_.detector_circles)
            r = r * RotMatrix(static_cast<float>(c.angle_deg * PI / 180.0), ImgCIFToInternal(c.axis));
        float rot1 = 0, rot2 = 0, rot3 = 0;
        PoniAnglesFromMatrix(r, 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));
    // The polarization degree is 2f - 1 of the fraction f the header states. Its sign is the plane:
    // positive where the plane contains the lab x axis (normal along y), negative where it contains y.
    // A laboratory tube says 0.5, unpolarized, where the plane does not matter.
    if (header0_.polarization_fraction.has_value()) {
        const Coord normal = ImgCIFToInternal(*header0_.polarization_normal);
        const float degree = 2.0f * static_cast<float>(*header0_.polarization_fraction) - 1.0f;
        dataset_->experiment.PolarizationFactor(std::fabs(normal.x) > std::fabs(normal.y) ? -degree : degree);
    }
    // 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("SMVReader").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());
    sweep::ForEachInOrder(files_.size(), [&](size_t i) {
        const auto h = smv::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};
    });

    // 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, "SMVReader");
    files_ = layout.files;
    dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name,
                                                   static_cast<float>(layout.start_deg),
                                                   static_cast<float>(layout.increment_deg),
                                                   header0_.spindle_axis.has_value()
                                                   ? ImgCIFToInternal(*header0_.spindle_axis)
                                                   : 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 JFJochSMVReader::GetNumberOfImages() const {
    return files_.size();
}

void JFJochSMVReader::Close() {
    files_.clear();
    dataset_.reset();
}

template <class Buffer>
CompressedImage JFJochSMVReader::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 = smv::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,
                              "SMV 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 JFJochSMVReader::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 JFJochSMVReader::HasImage(int64_t image_number) const {
    return image_number >= 0 && static_cast<size_t>(image_number) < files_.size()
           && !files_[image_number].empty();
}

bool JFJochSMVReader::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> JFJochSMVReader::ReadSpots(int64_t) const {
    return {};   // a raw SMV file stores no analysis results
}