diff --git a/apps/proxyncorr_gpu/main.cpp b/apps/proxyncorr_gpu/main.cpp index 20f1352..e0c3f11 100644 --- a/apps/proxyncorr_gpu/main.cpp +++ b/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(f) + 1); std::ofstream js((fs::path(output_dir) / (stem + ".json")).string()); diff --git a/include/cuncorr/dic_engine.h b/include/cuncorr/dic_engine.h index cb79f82..ce964be 100644 --- a/include/cuncorr/dic_engine.h +++ b/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: diff --git a/include/cuncorr/session.h b/include/cuncorr/session.h index abc0666..450ea89 100644 --- a/include/cuncorr/session.h +++ b/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; diff --git a/src/dic_engine.cpp b/src/dic_engine.cpp index 0001368..08a7f06 100644 --- a/src/dic_engine.cpp +++ b/src/dic_engine.cpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #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 seq_warp; // converged warp per grid point (last frame) + std::vector 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(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 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(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, std::greater> 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(gi) * g.width + gj], - static_cast(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(gi) * g.width + gj], + static_cast(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(gi) * g.width + gj], + static_cast(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(g.count()); + const double nan = std::numeric_limits::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(cr) * W + cc] == 0) return false; + return true; + }; + + std::vector> centers; + std::vector> icenters; + std::vector 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(center_r(gi)), + static_cast(center_c(gj))}); + icenters.push_back({center_r(gi), center_c(gj)}); + idxs.push_back(static_cast(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 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> cold_centers; + std::vector 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 diff --git a/src/session.cpp b/src/session.cpp index 3b7874b..8177714 100644 --- a/src/session.cpp +++ b/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