feat: sequence mode (cross-frame warm-start) + ROI/seed CLI

apps/proxyncorr_gpu/main.cpp

@@ -37,6 +37,10 @@ void print_usage(const char* argv0) {
         << "  --strain-radius N      Strain LS window radius (default 5)\n"
         << "  --threads N            Host orchestration threads (default 4)\n"
         << "  --ref PATH             Reference image (default: first frame)\n"
+        << "  --roi PATH             ROI mask image (grayscale > 0.5 = analysed)\n"
+        << "  --seed X,Y             Manual CPU seed at pixel (x=col, y=row); default auto\n"
+        << "  --seed-search N        Coarse seed search radius in px (default 15)\n"
+        << "  --sequence             Warm-start each frame from the previous (track growth)\n"
         << "  --no-strain            Skip strain computation\n"
         << "  --backend              Print the active compute backend and exit\n"
         << "  -h, --help             Show this help\n";
@@ -62,8 +66,8 @@ std::string frame_stem(int idx) {
 
 int main(int argc, char** argv) {
     cuncorr::SessionConfig cfg;
-    bool show_backend = false, do_strain = true;
-    std::string images_dir, output_dir, ref_path;
+    bool show_backend = false, do_strain = true, sequence = false;
+    std::string images_dir, output_dir, ref_path, roi_path;
 
     for (int i = 1; i < argc; ++i) {
         const std::string a = argv[i];
@@ -78,6 +82,16 @@ int main(int argc, char** argv) {
         else if (a == "--strain-radius") cfg.strain_radius = std::stoi(next("--strain-radius"));
         else if (a == "--threads") cfg.num_threads = std::stoi(next("--threads"));
         else if (a == "--ref") ref_path = next("--ref");
+        else if (a == "--roi") roi_path = next("--roi");
+        else if (a == "--seed-search") cfg.seed_search = std::stoi(next("--seed-search"));
+        else if (a == "--seed") {
+            const std::string s = next("--seed");
+            const auto comma = s.find(',');
+            if (comma == std::string::npos) { std::cerr << "--seed needs X,Y\n"; return 2; }
+            cfg.seed_x = std::stoi(s.substr(0, comma));
+            cfg.seed_y = std::stoi(s.substr(comma + 1));
+        }
+        else if (a == "--sequence") sequence = true;
         else if (a == "--no-strain") do_strain = false;
         else if (!a.empty() && a[0] == '-') { std::cerr << "unknown option: " << a << "\n"; return 2; }
         else if (images_dir.empty()) images_dir = a;
@@ -98,7 +112,14 @@ int main(int argc, char** argv) {
     else {
         frames.erase(std::remove(frames.begin(), frames.end(), ref_path), frames.end());
     }
+    // The ROI mask, if it lives inside the images dir, must not be treated as a frame.
+    if (!roi_path.empty())
+        frames.erase(std::remove(frames.begin(), frames.end(), roi_path), frames.end());
     std::cout << "reference: " << ref_path << "\ndeformed frames: " << frames.size() << "\n";
+    if (!roi_path.empty()) std::cout << "roi: " << roi_path << "\n";
+    if (cfg.seed_x >= 0) std::cout << "manual seed: (" << cfg.seed_x << "," << cfg.seed_y << ")\n";
+    else std::cout << "seed: auto\n";
+    if (sequence) std::cout << "mode: sequence (cross-frame warm-start)\n";
 
     const auto ref_img = cuncorr_cli::load_gray(ref_path);
     if (!ref_img.ok()) { std::cerr << "failed to load reference " << ref_path << "\n"; return 1; }
@@ -107,13 +128,22 @@ int main(int argc, char** argv) {
     cuncorr::NcorrSession session(cfg);
     session.set_reference(cuncorr::ImageBuffer(ref_img.data.data(), ref_img.width, ref_img.height, 1));
 
+    cuncorr_cli::GrayImage roi_img;
+    if (!roi_path.empty()) {
+        roi_img = cuncorr_cli::load_gray(roi_path);
+        if (!roi_img.ok()) { std::cerr << "failed to load roi " << roi_path << "\n"; return 1; }
+        session.set_roi(cuncorr::ImageBuffer(roi_img.data.data(), roi_img.width, roi_img.height, 1));
+    }
+
+    if (sequence) session.reset_sequence();
     int processed = 0, failed = 0;
     for (std::size_t f = 0; f < frames.size(); ++f) {
         const auto def = cuncorr_cli::load_gray(frames[f]);
         if (!def.ok()) { std::cerr << "  skip unreadable " << frames[f] << "\n"; ++failed; continue; }
 
-        const auto res = session.process_frame(
-            cuncorr::ImageBuffer(def.data.data(), def.width, def.height, 1));
+        const cuncorr::ImageBuffer def_buf(def.data.data(), def.width, def.height, 1);
+        const auto res = sequence ? session.process_frame_sequence(def_buf)
+                                  : session.process_frame(def_buf);
         const std::string stem = frame_stem(static_cast<int>(f) + 1);
 
         std::ofstream js((fs::path(output_dir) / (stem + ".json")).string());

include/cuncorr/dic_engine.h

@@ -53,8 +53,18 @@ class DicEngine {
     bool has_reference() const;
 
     /// Run DIC on @p def against the reference and return native Lagrangian displacements.
+    /// Stateless: every call seeds from scratch against the fixed reference.
     DICResult process(const ImageBuffer& def);
 
+    /// Sequence mode: like @ref process, but warm-starts each grid point from the previous
+    /// frame's converged warp (cross-frame seed propagation). This tracks a monotonically
+    /// growing displacement against the fixed reference without an ever-larger seed search,
+    /// mirroring CppNCorr's per-frame RGDIC. Call @ref reset_sequence before a new sequence.
+    DICResult process_sequence(const ImageBuffer& def);
+
+    /// Clear the cross-frame warm-start state (start a fresh sequence).
+    void reset_sequence();
+
     const BackendInfo& backend_info() const;
 
   private:

include/cuncorr/session.h

@@ -73,6 +73,17 @@ struct SessionConfig {
     int strain_radius = 5;     ///< Strain subregion radius in pixels.
     int num_threads = 4;       ///< Worker threads (CPU backend / host orchestration).
     bool debug = false;        ///< Verbose engine output.
+
+    // --- cuNCorr extensions (not part of the CppNCorr mirror) -------------------------
+    // Optional manual seed for the CPU reliability-guided path, in image pixel coords
+    // (seed_x = column, seed_y = row). -1 means auto-seed (scan for the first admissible,
+    // textured grid point). Ignored by the GPU batch path, which seeds every subset.
+    int seed_x = -1;           ///< Manual seed column (pixels); -1 = auto.
+    int seed_y = -1;           ///< Manual seed row (pixels); -1 = auto.
+    // Integer ZNCC coarse-search radius (px) used to lock a seed before IC-GN. Must cover
+    // the largest expected displacement (frames are matched against the fixed reference, so
+    // motion grows over a sequence). Larger = more robust but O(radius^2) per seed attempt.
+    int seed_search = 15;      ///< Coarse seed search radius in pixels.
 };
 
 /**
@@ -104,6 +115,14 @@ class NcorrSession {
     /// Push a deformed frame and run DIC. @throws std::logic_error if no reference set.
     DICResult process_frame(const ImageBuffer& def);
 
+    /// Sequence variant: warm-starts each grid point from the previous frame's result so a
+    /// growing displacement is tracked against the fixed reference without an ever-larger
+    /// seed search. Call @ref reset_sequence to begin a fresh sequence.
+    DICResult process_frame_sequence(const ImageBuffer& def);
+
+    /// Reset cross-frame warm-start state (start a new sequence).
+    void reset_sequence();
+
     /// @return true once a valid reference frame has been set.
     bool has_reference() const;
 

src/dic_engine.cpp

@@ -21,6 +21,7 @@
 #include <cmath>
 #include <limits>
 #include <queue>
+#include <stdexcept>
 #include <vector>
 
 #include "cuncorr/icgn.h"
@@ -79,10 +80,14 @@ struct DicEngine::State {
     bool has_ref = false;
 
     // engine tuning (cuNCorr internals; not part of the session.h contract)
-    int seed_search = 15;         // integer ZNCC search radius for the seed only
     double accept_cznssd = 1.0;   // residual gate for propagation acceptance
     IcgnConfig icgn;              // bicubic, 50 iters, 1e-4
 
+    // --- cross-frame warm-start state for process_sequence() --------------------------
+    bool seq_started = false;
+    std::vector<WarpParams> seq_warp;  // converged warp per grid point (last frame)
+    std::vector<char> seq_valid;       // 1 where seq_warp holds a usable warm-start
+
     explicit State(const SessionConfig& c)
         : cfg(c), backend(make_default_backend()), backend_info(backend->info()) {}
 };
@@ -162,7 +167,7 @@ DICResult DicEngine::process(const ImageBuffer& def) {
                     idxs.push_back(static_cast<std::size_t>(gi) * g.width + gj);
                 }
         if (!centers.empty()) {
-            const auto seeds = coarse_batch_cuda(st_->ref, cur, subset, icenters, st_->seed_search);
+            const auto seeds = coarse_batch_cuda(st_->ref, cur, subset, icenters, st_->cfg.seed_search);
             std::vector<WarpParams> inits(centers.size());
             for (std::size_t i = 0; i < inits.size(); ++i)
                 inits[i] = seeds[i].valid ? seeds[i].p : WarpParams{};
@@ -190,7 +195,7 @@ DICResult DicEngine::process(const ImageBuffer& def) {
     auto try_seed = [&](int gi, int gj) -> bool {
         const std::size_t idx = static_cast<std::size_t>(gi) * g.width + gj;
         const CoarseGuess cg = coarse_init(st_->ref, cur, subset, center_r(gi), center_c(gj),
-                                           st_->seed_search);
+                                           st_->cfg.seed_search);
         const WarpParams init = cg.valid ? cg.p : WarpParams{};
         const IcgnResult res = opt.optimize(center_r(gi), center_c(gj), init);
         if (res.valid && res.converged && res.cznssd < st_->accept_cznssd) {
@@ -209,14 +214,32 @@ DICResult DicEngine::process(const ImageBuffer& def) {
     std::priority_queue<QItem, std::vector<QItem>, std::greater<QItem>> queue;
 
     bool seeded = false;
-    for (int gi = 0; gi < g.height && !seeded; ++gi)
-        for (int gj = 0; gj < g.width && !seeded; ++gj)
-            if (admissible(gi, gj) && try_seed(gi, gj)) {
-                queue.emplace(r.corrcoef[static_cast<std::size_t>(gi) * g.width + gj],
-                              static_cast<std::size_t>(gi) * g.width + gj);
-                seeded = true;
-            }
-    if (!seeded) { r.message = "no valid seed found (check ROI / texture / radius)"; return r; }
+    if (st_->cfg.seed_x >= 0 && st_->cfg.seed_y >= 0) {
+        // Manual seed: snap the requested pixel (x=col, y=row) to the nearest grid point.
+        const int gi = (st_->cfg.seed_y + sf / 2) / sf;
+        const int gj = (st_->cfg.seed_x + sf / 2) / sf;
+        if (gi >= 0 && gi < g.height && gj >= 0 && gj < g.width && admissible(gi, gj) &&
+            try_seed(gi, gj)) {
+            queue.emplace(r.corrcoef[static_cast<std::size_t>(gi) * g.width + gj],
+                          static_cast<std::size_t>(gi) * g.width + gj);
+            seeded = true;
+        }
+        if (!seeded) {
+            r.message = "manual seed (" + std::to_string(st_->cfg.seed_x) + "," +
+                        std::to_string(st_->cfg.seed_y) + ") failed (off-grid / outside ROI / "
+                        "untextured / no convergence)";
+            return r;
+        }
+    } else {
+        for (int gi = 0; gi < g.height && !seeded; ++gi)
+            for (int gj = 0; gj < g.width && !seeded; ++gj)
+                if (admissible(gi, gj) && try_seed(gi, gj)) {
+                    queue.emplace(r.corrcoef[static_cast<std::size_t>(gi) * g.width + gj],
+                                  static_cast<std::size_t>(gi) * g.width + gj);
+                    seeded = true;
+                }
+        if (!seeded) { r.message = "no valid seed found (check ROI / texture / radius)"; return r; }
+    }
 
     // --- reliability-guided flood fill -------------------------------------------------
     const int dgi[4] = {-1, 1, 0, 0};
@@ -258,4 +281,125 @@ DICResult DicEngine::process(const ImageBuffer& def) {
     return r;
 }
 
+void DicEngine::reset_sequence() {
+    st_->seq_started = false;
+    st_->seq_warp.clear();
+    st_->seq_valid.clear();
+}
+
+DICResult DicEngine::process_sequence(const ImageBuffer& def) {
+    if (!st_->has_ref) throw std::logic_error("DicEngine::process_sequence: no reference set");
+
+    DICResult r;
+    if (!def.valid() || def.width != st_->ref.w || def.height != st_->ref.h) {
+        r.message = "deformed frame geometry does not match reference";
+        return r;
+    }
+
+    const Image2Df cur = to_image(def);
+    const int sf = st_->cfg.scalefactor > 0 ? st_->cfg.scalefactor : 1;
+    const int radius = st_->cfg.subregion_radius;
+    const SubsetOffsets subset = make_subset(radius, SubsetShape::CIRCLE);
+    const int margin = radius + 2;
+    const int H = st_->ref.h, W = st_->ref.w;
+
+    const GridGeometry g = reduced_grid(W, H, st_->cfg);
+    r.width = g.width;
+    r.height = g.height;
+    const std::size_t N = static_cast<std::size_t>(g.count());
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+    r.u.assign(N, nan);
+    r.v.assign(N, nan);
+    r.corrcoef.assign(N, nan);
+    if (N == 0) { r.message = "empty grid"; return r; }
+
+    if (!st_->seq_started || st_->seq_warp.size() != N) {
+        st_->seq_warp.assign(N, WarpParams{});
+        st_->seq_valid.assign(N, 0);
+        st_->seq_started = true;
+    }
+
+    auto center_r = [&](int gi) { return gi * sf; };
+    auto center_c = [&](int gj) { return gj * sf; };
+    auto admissible = [&](int gi, int gj) {
+        const int cr = center_r(gi), cc = center_c(gj);
+        if (cr < margin || cr > H - 1 - margin || cc < margin || cc > W - 1 - margin) return false;
+        if (!st_->roi.empty() && st_->roi[static_cast<std::size_t>(cr) * W + cc] == 0) return false;
+        return true;
+    };
+
+    std::vector<std::array<double, 2>> centers;
+    std::vector<std::array<int, 2>> icenters;
+    std::vector<std::size_t> idxs;
+    for (int gi = 0; gi < g.height; ++gi)
+        for (int gj = 0; gj < g.width; ++gj)
+            if (admissible(gi, gj)) {
+                centers.push_back({static_cast<double>(center_r(gi)),
+                                   static_cast<double>(center_c(gj))});
+                icenters.push_back({center_r(gi), center_c(gj)});
+                idxs.push_back(static_cast<std::size_t>(gi) * g.width + gj);
+            }
+    if (centers.empty()) { r.message = "no admissible points"; return r; }
+
+    // Warm-start init per point: previous frame's converged warp where available, else a
+    // fresh coarse seed (needed on frame 0 and to recover points that dropped out).
+    std::vector<WarpParams> inits(centers.size());
+
+#ifdef CUNCORR_WITH_CUDA
+    if (st_->backend_info.kind == Backend::CUDA && st_->backend_info.device_available) {
+        // Coarse-seed only the points without a warm-start (frame 0 -> all of them).
+        std::vector<std::array<int, 2>> cold_centers;
+        std::vector<std::size_t> cold_pos;
+        for (std::size_t i = 0; i < idxs.size(); ++i) {
+            if (st_->seq_valid[idxs[i]]) inits[i] = st_->seq_warp[idxs[i]];
+            else { cold_centers.push_back(icenters[i]); cold_pos.push_back(i); }
+        }
+        if (!cold_centers.empty()) {
+            const auto seeds = coarse_batch_cuda(st_->ref, cur, subset, cold_centers,
+                                                 st_->cfg.seed_search);
+            for (std::size_t k = 0; k < cold_pos.size(); ++k)
+                inits[cold_pos[k]] = seeds[k].valid ? seeds[k].p : WarpParams{};
+        }
+        const auto results = icgn_batch_cuda(st_->ref, cur, subset, st_->icgn, centers, inits);
+        for (std::size_t i = 0; i < results.size(); ++i) {
+            const auto& rr = results[i];
+            const std::size_t idx = idxs[i];
+            if (rr.valid && rr.converged && rr.cznssd < st_->accept_cznssd) {
+                r.u[idx] = rr.p.u; r.v[idx] = rr.p.v; r.corrcoef[idx] = rr.cznssd;
+                st_->seq_warp[idx] = rr.p; st_->seq_valid[idx] = 1;
+            } else {
+                st_->seq_valid[idx] = 0;  // dropped out; re-seed coarsely next frame
+            }
+        }
+        r.valid = true;
+        r.message = "ok (gpu sequence, backend=" + st_->backend_info.name + ")";
+        return r;
+    }
+#endif
+
+    // CPU reference path: warm-start each point, coarse-seed the cold ones.
+    IcgnOptimizer opt(st_->ref, cur, subset, st_->icgn);
+    for (std::size_t i = 0; i < idxs.size(); ++i) {
+        const std::size_t idx = idxs[i];
+        WarpParams init;
+        if (st_->seq_valid[idx]) {
+            init = st_->seq_warp[idx];
+        } else {
+            const CoarseGuess cg = coarse_init(st_->ref, cur, subset, icenters[i][0],
+                                               icenters[i][1], st_->cfg.seed_search);
+            init = cg.valid ? cg.p : WarpParams{};
+        }
+        const IcgnResult res = opt.optimize(icenters[i][0], icenters[i][1], init);
+        if (res.valid && res.converged && res.cznssd < st_->accept_cznssd) {
+            r.u[idx] = res.p.u; r.v[idx] = res.p.v; r.corrcoef[idx] = res.cznssd;
+            st_->seq_warp[idx] = res.p; st_->seq_valid[idx] = 1;
+        } else {
+            st_->seq_valid[idx] = 0;
+        }
+    }
+    r.valid = true;
+    r.message = "ok (cpu sequence, backend=" + st_->backend_info.name + ")";
+    return r;
+}
+
 }  // namespace cuncorr

src/session.cpp

@@ -25,6 +25,10 @@ NcorrSession& NcorrSession::operator=(NcorrSession&&) noexcept = default;
 void NcorrSession::set_reference(const ImageBuffer& ref) { impl_->engine.set_reference(ref); }
 void NcorrSession::set_roi(const ImageBuffer& roi_mask) { impl_->engine.set_roi(roi_mask); }
 DICResult NcorrSession::process_frame(const ImageBuffer& def) { return impl_->engine.process(def); }
+DICResult NcorrSession::process_frame_sequence(const ImageBuffer& def) {
+    return impl_->engine.process_sequence(def);
+}
+void NcorrSession::reset_sequence() { impl_->engine.reset_sequence(); }
 bool NcorrSession::has_reference() const { return impl_->engine.has_reference(); }
 
 }  // namespace cuncorr