Files
Jungfraujoch/reader/HDF5ImageSource.cpp
T
leonarski_fandClaude Opus 5 1d16dcd2b9
Build Packages / build:windows:nocuda (push) Successful in 16m8s
Build Packages / build:windows:cuda (push) Successful in 18m58s
Build Packages / build:viewer-tgz:cpu (push) Successful in 20m35s
Build Packages / build:viewer-tgz:cuda (push) Successful in 22m31s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 25m9s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 25m6s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 28m57s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 28m58s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 28m58s
Build Packages / XDS test (durin plugin) (push) Successful in 12m3s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 22m24s
Build Packages / build:rpm (rocky9) (push) Successful in 21m45s
Build Packages / Generate python client (push) Successful in 53s
Build Packages / build:rpm (rocky8) (push) Successful in 26m9s
Build Packages / Create release (push) Skipped
Build Packages / Build documentation (push) Successful in 1m37s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 25m34s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 22m0s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 10m53s
Build Packages / XDS test (neggia plugin) (push) Successful in 9m29s
Build Packages / DIALS test (push) Successful in 23m40s
Build Packages / Unit tests (push) Successful in 1h20m1s
Read a chunk without zeroing it first, and hold no frame the fused decoder never writes
Three costs before and around the image loop.

Every image allocated a fresh buffer for its compressed chunk and resized it, which
value-initialises, and the read then overwrote every byte. At a few megabytes a chunk
the allocation is large enough to be mapped rather than reused, so the zeroing was
page-fault bound and cost more than the read it preceded - twenty gigabytes of it
over a long sweep. The buffer now uses an allocator that does not construct, and the
two HDF5 read paths are templated on the allocator so every existing caller compiles
unchanged. The rebind is deliberate: without it the vector base rebinds to the
default allocator and the zeroing quietly returns.

The bitshuffle decoder allocated a whole uncompressed frame in its constructor -
seventy megabytes a worker, five hundred and fifty across the loop - for the route
that decodes the shuffled image separately. That route is taken only when a
bitshuffle block is too large for the fused kernel, which neither writer this
pipeline reads produces, so on a real frame the buffer is allocated, never touched,
and freed. It is now allocated where it is used. The comment two lines below already
warned against sizing a buffer from the uncompressed size; the line above it had not
been given the same treatment.

The first call into cuFFT pays the library's one-time initialisation, and it landed
in the middle of the first pass with nothing to overlap it. It is now forced on a
background thread at startup, alongside the file open and the mapping build, in the
manner the shadow finder already uses.

Finally, the detector mask was copied into the start message whether or not a file
would carry it, which a merging run does not. It is filled where a writer is
constructed - both places one is constructed, the second being the fallback that
writes a process file when nothing indexed.

Faster on eleven of thirty-eight crystals and slower on none; the whole rotation test
set falls from four minutes thirty to four minutes seventeen, with each binary
repeating itself to within half a per cent. Space groups thirty-five of thirty-eight
and no failures throughout, and every column of the comparison table is identical.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EGpGdgmJ8MyY9pCGWjktyi
2026-08-25 01:08:38 +02:00

160 lines
6.0 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "HDF5ImageSource.h"
#include "../common/JFJochException.h"
#ifdef _WIN32
#include <windows.h>
#else
#include <fcntl.h>
#include <unistd.h>
#endif
// Positional reads: pread() on POSIX, ReadFile() with an OVERLAPPED offset on Windows. Both take the
// offset as an argument instead of moving a shared file position, which is what lets every worker
// thread read through one handle at the same time.
HDF5ImageSource::RawFile::RawFile(const std::string &path) {
#ifdef _WIN32
HANDLE h = CreateFileA(path.c_str(), GENERIC_READ, FILE_SHARE_READ | FILE_SHARE_WRITE, nullptr,
OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, nullptr);
handle_ = (h == INVALID_HANDLE_VALUE) ? -1 : reinterpret_cast<intptr_t>(h);
#else
handle_ = ::open(path.c_str(), O_RDONLY);
#endif
}
HDF5ImageSource::RawFile::~RawFile() {
if (handle_ == -1)
return;
#ifdef _WIN32
CloseHandle(reinterpret_cast<HANDLE>(handle_));
#else
::close(static_cast<int>(handle_));
#endif
}
void HDF5ImageSource::RawFile::ReadAt(void *dst, size_t size, uint64_t address) const {
auto *out = static_cast<uint8_t *>(dst);
size_t done = 0;
while (done < size) {
#ifdef _WIN32
OVERLAPPED ov{};
ov.Offset = static_cast<DWORD>((address + done) & 0xFFFFFFFFULL);
ov.OffsetHigh = static_cast<DWORD>((address + done) >> 32);
DWORD got = 0;
const bool ok = ReadFile(reinterpret_cast<HANDLE>(handle_), out + done,
static_cast<DWORD>(size - done), &got, &ov);
const long long n = ok ? static_cast<long long>(got) : -1;
#else
const long long n = ::pread(static_cast<int>(handle_), out + done, size - done, address + done);
#endif
if (n <= 0)
throw JFJochException(JFJochExceptionCategory::HDF5, "Error reading image chunk from file");
done += static_cast<size_t>(n);
}
}
void HDF5ImageSource::Configure(HDF5ImageLocator::Layout layout) {
dataset_cache_.clear();
locator_.Configure(std::move(layout));
}
void HDF5ImageSource::Clear() {
dataset_cache_.clear();
locator_.Clear();
}
HDF5ImageLocator::Location HDF5ImageSource::Resolve(int64_t global) const {
return locator_.Resolve(global);
}
StoredPixelFormat HDF5ImageSource::GetStoredPixelFormat() const {
auto loc = locator_.Resolve(0);
HDF5DataSet dataset(*loc.file, "/entry/data/data");
HDF5DataType datatype(dataset);
return {static_cast<int64_t>(datatype.GetElemSize()) * 8, datatype.IsSigned()};
}
std::vector<HDF5DataSourceMessage> HDF5ImageSource::GetSourceMapping(uint64_t first_image,
std::optional<uint64_t> image_count,
uint64_t total_images,
uint64_t stride) const {
return locator_.GetSourceMapping(first_image, image_count, total_images, stride);
}
const HDF5ImageSource::OpenDataset &
HDF5ImageSource::GetDataset(const HDF5ImageLocator::Location &loc) const {
if (auto it = dataset_cache_.find(loc.file.get()); it != dataset_cache_.end())
return it->second;
OpenDataset entry;
entry.file = loc.file;
entry.dataset = std::make_unique<HDF5DataSet>(*loc.file, "/entry/data/data");
HDF5DataSpace dataspace(*entry.dataset);
HDF5DataType datatype(*entry.dataset);
HDF5Dcpl dcpl(*entry.dataset);
if (dataspace.GetNumOfDimensions() != 3)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"/entry/data/data dataset must be 3D");
if (datatype.IsFloat())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Float datasets not supported at this time");
auto dim = dataspace.GetDimensions();
entry.height = dim[1];
entry.width = dim[2];
entry.mode = CalcImageMode(datatype.GetElemSize(), datatype.IsFloat(), datatype.IsSigned());
auto chunk_size = dcpl.GetChunking();
entry.direct_chunk = (chunk_size.size() == 3) && (chunk_size[0] == 1)
&& (chunk_size[1] == dim[1]) && (chunk_size[2] == dim[2]);
if (entry.direct_chunk)
entry.algorithm = dcpl.GetCompression();
if (entry.direct_chunk && !loc.path.empty()) {
entry.raw = std::make_shared<RawFile>(loc.path);
if (!entry.raw->IsOpen())
entry.raw.reset();
hid_t fcpl = H5Fget_create_plist(loc.file->GetID());
if (fcpl >= 0) {
hsize_t user_block = 0;
if (H5Pget_userblock(fcpl, &user_block) >= 0)
entry.user_block = user_block;
H5Pclose(fcpl);
}
}
return dataset_cache_.emplace(loc.file.get(), std::move(entry)).first->second;
}
std::optional<HDF5ImageSource::DirectChunk>
HDF5ImageSource::PrepareDirectRead(const HDF5ImageLocator::Location &loc) const {
const auto &ds = GetDataset(loc);
if (!ds.raw)
return {};
const hsize_t coord[3] = {static_cast<hsize_t>(loc.local_index), 0, 0};
unsigned filter_mask = 0;
haddr_t address = HADDR_UNDEF;
hsize_t size = 0;
if (H5Dget_chunk_info_by_coord(ds.dataset->GetID(), coord, &filter_mask, &address, &size) < 0)
return {};
// A chunk nobody ever wrote has no address and no bytes; only HDF5 knows it reads as the fill
// value, so hand those back to it.
if (address == HADDR_UNDEF || size == 0)
return {};
return DirectChunk{ds.raw, ds.user_block + address, static_cast<uint32_t>(size),
ds.width, ds.height, ds.mode, ds.algorithm};
}
CompressedImage HDF5ImageSource::ReadDirect(RawByteBuffer &buffer, const DirectChunk &chunk) {
buffer.resize(chunk.size);
chunk.file->ReadAt(buffer.data(), chunk.size, chunk.address);
return {buffer.data(), buffer.size(), chunk.width, chunk.height, chunk.mode, chunk.algorithm};
}