Move the atomic-model code to image_analysis/structure_refinement/ and WriteModel to writer/

A pure move. ModelValidation, RigidBodyRefine, RigidBodyGPU, ModelFFT, ModelGrid, ModelScaling, ModelMaskGPU, ModelScaleGPU and SigmaA - everything that works on an atomic model - become the JFJochStructureRefinement library, linked by JFJochImageAnalysis. WriteModel (the placed-model mmCIF/PDB writer) goes to writer/ as its own small JFJochModelWriter target, so JFJochWriter, which a writer-only build compiles, does not gain a gemmi dependency. Only include paths and CMake lists change. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
2026-10-07 14:05:37 +02:00
parent 08ca1d5764
commit 9ad92b6bfe
34 changed files with 73 additions and 69 deletions
@@ -0,0 +1,112 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ModelFFT.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <complex>
+#include <map>
+#include <mutex>
+
+#include <fftw3.h>
+
+#include "gemmi/fail.hpp"
+
+#include "../../common/FFTWPlannerLock.h"
+
+namespace {
+
+// The r2c plan (or, inverse, the c2r plan) for an (nu, nv, nw) map stored u fastest, halving w - the
+// layout gemmi's own transforms read and write. FFTW's planner is not thread-safe, so plans are made
+// under the process-wide planner lock (which also guards this cache); executing one on new arrays
+// (fftwf_execute_dft_r2c / _c2r) is, which is how every caller uses it.
+fftwf_plan PlanFor(int nu, int nv, int nw, bool inverse) {
+    static std::map<std::array<int, 4>, fftwf_plan> plans;
+    std::lock_guard lock(FFTWPlannerMutex());
+    const std::array<int, 4> key{nu, nv, nw, inverse};
+    if (const auto it = plans.find(key); it != plans.end())
+        return it->second;
+    float *real = fftwf_alloc_real(static_cast<size_t>(nu) * nv * nw);
+    fftwf_complex *cplx = fftwf_alloc_complex(static_cast<size_t>(nu) * nv * (nw / 2 + 1));
+    // The last dimension is the one r2c halves; the strides put u fastest on both sides.
+    fftwf_iodim dims[3] = {{nu, 1, 1}, {nv, nu, nu}, {nw, nu * nv, nu * nv}};
+    fftwf_plan plan = nullptr;
+    if (real != nullptr && cplx != nullptr)
+        plan = inverse ? fftwf_plan_guru_dft_c2r(3, dims, 0, nullptr, cplx, real, FFTW_ESTIMATE)
+                       : fftwf_plan_guru_dft_r2c(3, dims, 0, nullptr, real, cplx, FFTW_ESTIMATE);
+    fftwf_free(real);
+    fftwf_free(cplx);
+    if (plan == nullptr)
+        gemmi::fail("ModelFFT: no FFTW plan for the grid");
+    plans.emplace(key, plan);
+    return plan;
+}
+
+}  // namespace
+
+gemmi::FPhiGrid<float> MapToFPhi(const gemmi::Grid<float> &map) {
+    if (map.axis_order == gemmi::AxisOrder::ZYX)
+        gemmi::fail("MapToFPhi(): ZYX order is not supported");
+    gemmi::FPhiGrid<float> hkl;
+    hkl.unit_cell = map.unit_cell;
+    hkl.spacegroup = map.spacegroup;
+    hkl.axis_order = map.axis_order;
+    hkl.half_l = true;
+    hkl.set_size_without_checking(map.nu, map.nv, map.nw / 2 + 1);
+    const float norm = static_cast<float>(map.unit_cell.volume / map.point_count());
+
+    const fftwf_plan plan = PlanFor(map.nu, map.nv, map.nw, false);
+    // Buffers from fftwf_malloc: a plan may only be executed on arrays aligned as the ones it was made
+    // on were, which the vectors' own allocations do not promise.
+    float *in = fftwf_alloc_real(map.data.size());
+    fftwf_complex *out = fftwf_alloc_complex(hkl.data.size());
+    if (in == nullptr || out == nullptr) {
+        fftwf_free(in);
+        fftwf_free(out);
+        gemmi::fail("MapToFPhi(): cannot allocate the FFT buffers");
+    }
+    std::copy(map.data.begin(), map.data.end(), in);
+    fftwf_execute_dft_r2c(plan, in, out);
+    // Conjugated, as gemmi does: FFTW's forward transform carries exp(-2 pi i h.x), and a structure
+    // factor is the sum over exp(+2 pi i h.x).
+    for (size_t i = 0; i < hkl.data.size(); ++i)
+        hkl.data[i] = std::complex<float>(out[i][0] * norm, -out[i][1] * norm);
+    fftwf_free(in);
+    fftwf_free(out);
+    return hkl;
+}
+
+gemmi::Grid<float> MapFromFPhi(const gemmi::FPhiGrid<float> &hkl) {
+    if (hkl.axis_order == gemmi::AxisOrder::ZYX || !hkl.half_l)
+        gemmi::fail("MapFromFPhi(): only an XYZ half-l grid is supported");
+    gemmi::Grid<float> map;
+    map.spacegroup = hkl.spacegroup;
+    map.unit_cell = hkl.unit_cell;
+    map.set_size(hkl.nu, hkl.nv, 2 * (hkl.nw - 1));
+    map.axis_order = hkl.axis_order;
+    const float norm = static_cast<float>(1.0 / hkl.unit_cell.volume);
+
+    const fftwf_plan plan = PlanFor(map.nu, map.nv, map.nw, true);
+    fftwf_complex *in = fftwf_alloc_complex(hkl.data.size());
+    float *out = fftwf_alloc_real(map.data.size());
+    if (in == nullptr || out == nullptr) {
+        fftwf_free(in);
+        fftwf_free(out);
+        gemmi::fail("MapFromFPhi(): cannot allocate the FFT buffers");
+    }
+    // Conjugated, as gemmi does: rho(x) = 1/V sum F exp(-2 pi i h.x), and FFTW's backward transform
+    // carries exp(+2 pi i h.x). A missing coefficient (NaN) counts as zero, also as gemmi does.
+    for (size_t i = 0; i < hkl.data.size(); ++i) {
+        const std::complex<float> x = hkl.data[i];
+        in[i][0] = std::isnan(x.imag()) ? 0.0f : x.real();
+        in[i][1] = std::isnan(x.imag()) ? 0.0f : -x.imag();
+    }
+    fftwf_execute_dft_c2r(plan, in, out);
+    for (size_t i = 0; i < map.data.size(); ++i)
+        map.data[i] = out[i] * norm;
+    fftwf_free(in);
+    fftwf_free(out);
+    return map;
+}