ncorr.h

Go to the documentation of this file.

/* 
 * File:   ncorr.h
 * Author: justin
 *
 * Created on May 12, 2015, 1:33 AM
 */
 
#ifndef NCORR_H
#define NCORR_H
 
#include "Array2D.h"
#include "Image2D.h"
#include "ROI2D.h"
#include "Data2D.h"
#include "Disp2D.h"
#include "Strain2D.h"
 
namespace ncorr {
 
namespace details {
    // Nonlinear optimization ------------------------------------------------//
    class nloptimizer_base {         
        public:      
            typedef std::ptrdiff_t                              difference_type;    
            typedef std::pair<difference_type,difference_type>           coords;  
                        
            // Rule of 5 and destructor --------------------------------------//
            nloptimizer_base() noexcept = default;
            nloptimizer_base(const nloptimizer_base&) = default;
            nloptimizer_base(nloptimizer_base&&) = default;
            nloptimizer_base& operator=(const nloptimizer_base&) = default;  
            nloptimizer_base& operator=(nloptimizer_base&&) = default;
            virtual ~nloptimizer_base() noexcept = default;
            
            // Additional Constructors ---------------------------------------//         
            nloptimizer_base(difference_type order, difference_type num_params) : 
                grad_buf(order,1), hess_buf(order,order), params(num_params,1) { }
                         
            // Arithmetic operations -----------------------------------------//
            std::pair<const Array2D<double>&, bool> global(const Array2D<double>&) const;
            std::pair<const Array2D<double>&, bool> operator()(const Array2D<double>&) const;           
            
        protected:       
            // Arithmetic operations -----------------------------------------//
            virtual bool initial_guess() const = 0;
            virtual bool iterative_search() const = 0;
            virtual bool newton() const = 0;
                         
            // Utility -------------------------------------------------------//
            void chk_input_params_size(const Array2D<double>&) const;
            
            mutable Array2D<double> grad_buf; 
            mutable Array2D<double> hess_buf; 
            mutable Array2D<double> params; 
            double cutoff_norm = 1e-6;
            difference_type cutoff_iterations = 100;
    };
    
    class disp_nloptimizer final : public nloptimizer_base {
        public:      
            typedef nloptimizer_base::difference_type           difference_type;   
            typedef nloptimizer_base::coords                             coords;
                        
            // Rule of 5 and destructor --------------------------------------//
            disp_nloptimizer() noexcept : region_idx() { }
            disp_nloptimizer(const disp_nloptimizer&) = default;
            disp_nloptimizer(disp_nloptimizer&&) = default;
            disp_nloptimizer& operator=(const disp_nloptimizer&) = default;  
            disp_nloptimizer& operator=(disp_nloptimizer&&) = default;
            ~disp_nloptimizer() noexcept = default;
            
            // Additional Constructors ---------------------------------------//         
            // Note: params = {p1_new, p2_new, p1_old, p2_old, v_old, u_old, dv_dp1_old, dv_dp2_old, du_dp1_old, du_dp2_old, dist, grad_norm}
            disp_nloptimizer(const Disp2D &disp, difference_type region_idx, INTERP interp_type) : 
                nloptimizer_base(2, 12), disp(disp), region_idx(region_idx), disp_interp(disp.get_nlinfo_interpolator(region_idx, interp_type)) { }
                         
        private:       
            // Arithmetic operations -----------------------------------------//
            bool initial_guess() const override;
            bool iterative_search() const override;
            bool newton() const override;
            
            Disp2D disp;                             // immutable - Disp2D has pointer semantics
            difference_type region_idx;
            Disp2D::nlinfo_interpolator disp_interp; // Have copy
    };
    
    class subregion_nloptimizer final : public nloptimizer_base {
        public:      
            typedef nloptimizer_base::difference_type           difference_type;   
            typedef nloptimizer_base::coords                             coords;
                        
            // Rule of 5 and destructor --------------------------------------//
            subregion_nloptimizer() noexcept : scalefactor(), ref_template_avg(), ref_template_ssd_inv(), last_iteration_count(0) { }
            subregion_nloptimizer(const subregion_nloptimizer&) = default;
            subregion_nloptimizer(subregion_nloptimizer&&) = default;
            subregion_nloptimizer& operator=(const subregion_nloptimizer&) = default;  
            subregion_nloptimizer& operator=(subregion_nloptimizer&&) = default;
            ~subregion_nloptimizer() noexcept = default;
 
            // getters -------------------------------------------------------//
            const ROI2D::contig_subregion_generator& get_subregion_gen() const { return this->subregion_gen; }
            int get_last_iteration_count() const { return last_iteration_count; }
            
            
            // Additional Constructors ---------------------------------------//         
            // Note: params = {p1, p2, v, u, dv_dp1, dv_dp2, du_dp1, du_dp2, corr_coef, diff_norm}
            subregion_nloptimizer(const Array2D<double>&, const Array2D<double>&, const ROI2D&, difference_type, INTERP, SUBREGION, difference_type);
                         
        private:      
            // Arithmetic operations -----------------------------------------//
            bool initial_guess() const override;
            bool iterative_search() const override;
            bool newton() const override;
            
            std::shared_ptr<Array2D<double>> A_ref_ptr;                // Allows R-value arrays; immutable
            std::shared_ptr<Array2D<double>> A_cur_ptr;                // Allows R-value arrays; immutable
            difference_type scalefactor;
            Array2D<double>::interpolator A_cur_interp;                // Have copy
            ROI2D::contig_subregion_generator subregion_gen;           // Have copy
            // Buffers for NCC:
            std::shared_ptr<Array2D<double>> A_cur_cumsum_p1_ptr;      // immutable
            std::shared_ptr<Array2D<double>> A_cur_pow_cumsum_p1_ptr;  // immutable
            mutable Array2D<double> A_ref_template;                    // Have copy 
            // Buffers for inverse compositional gauss newton method:
            std::shared_ptr<Array2D<double>> A_dref_dp1_ptr;           // immutable
            std::shared_ptr<Array2D<double>> A_dref_dp2_ptr;           // immutable
            mutable double ref_template_avg;
            mutable double ref_template_ssd_inv;
            // Steepest descent images:
            mutable Array2D<double> A_dref_dv;                         // Have copy 
            mutable Array2D<double> A_dref_du;                         // Have copy 
            mutable Array2D<double> A_dref_dv_dp1;                     // Have copy 
            mutable Array2D<double> A_dref_dv_dp2;                     // Have copy 
            mutable Array2D<double> A_dref_du_dp1;                     // Have copy 
            mutable Array2D<double> A_dref_du_dp2;                     // Have copy 
            // Linsolver for hessian
            mutable Array2D<double>::linsolver hess_linsolver;         // Have copy 
            // Cur template buffer
            mutable Array2D<double> A_cur_template;                    // Have copy 
            // Iteration tracking for MATLAB-style saturation check
            mutable int last_iteration_count;                          // Track actual iterations
    };        
}
 
// Interface functions -------------------------------------------------------//
// ROI update modes:
// - SKIP_ALL: Original behavior - if ANY boundary point returns NaN, the entire
//             boundary is discarded (returns empty ROI on failure)
// - SKIP_INVALID: Skip individual NaN points and out-of-bounds points, keeping
//                 valid points. Only returns empty if ALL points are invalid.
enum class ROI_UPDATE_MODE { SKIP_ALL, SKIP_INVALID };
 
// Accumulation modes:
// - ON_THE_FLY: Original C++ behavior - accumulate displacements during analysis.
//               Fast but can fail when ROI changes shape during updates.
// - POST_PROCESS: MATLAB-like approach - store individual step displacements,
//                 then add them all at the end with their matching ROIs.
//                 More robust when ROI updates occur frequently.
enum class ACCUMULATION_MODE { ON_THE_FLY, POST_PROCESS };
 
ROI2D update(const ROI2D&, const Disp2D&, INTERP, ROI_UPDATE_MODE mode = ROI_UPDATE_MODE::SKIP_ALL);
ROI2D matlab_update_roi(const ROI2D&, const Disp2D&, INTERP, ROI2D::difference_type radius);
 
Data2D update(const Data2D&, const Disp2D&, INTERP, ROI_UPDATE_MODE mode = ROI_UPDATE_MODE::SKIP_ALL);
 
// Original add function - combines displacements using first disp's ROI
Disp2D add(const std::vector<Disp2D>&, INTERP);
 
// MATLAB-style add function - combines displacements using each disp's own ROI for bounds checking
Disp2D add_with_rois(const std::vector<Disp2D>& disps, const std::vector<ROI2D>& rois, INTERP interp_type);
 
// DIC_analysis --------------------------------------------------------------//
Disp2D RGDIC(const Array2D<double>&, const Array2D<double>&, const ROI2D&, ROI2D::difference_type, INTERP, SUBREGION, ROI2D::difference_type, ROI2D::difference_type, double, bool);
 
enum class DIC_analysis_config { NO_UPDATE, KEEP_MOST_POINTS, REMOVE_BAD_POINTS };
 
struct DIC_analysis_input final {
    typedef ROI2D::difference_type                              difference_type;
        
    // Rule of 5 and destructor ----------------------------------------------//    
    DIC_analysis_input() : scalefactor(), interp_type(), subregion_type(), r(), num_threads(), cutoff_corrcoef(), update_corrcoef(), prctile_corrcoef(), roi_update_mode(ROI_UPDATE_MODE::SKIP_ALL), accumulation_mode(ACCUMULATION_MODE::ON_THE_FLY), save_disps_steps(false), debug() { }
    DIC_analysis_input(const DIC_analysis_input&) = default;
    DIC_analysis_input(DIC_analysis_input&&) = default;
    DIC_analysis_input& operator=(const DIC_analysis_input&) = default;
    DIC_analysis_input& operator=(DIC_analysis_input&&) = default;
    ~DIC_analysis_input() noexcept = default;
    
    // Additional constructors -----------------------------------------------//    
    DIC_analysis_input(const std::vector<Image2D> &imgs,
                       const ROI2D &roi,
                       difference_type scalefactor, 
                       INTERP interp_type,
                       SUBREGION subregion_type,
                       difference_type r,
                       difference_type num_threads,
                       double cutoff_corrcoef,
                       double update_corrcoef,
                       double prctile_corrcoef,
                       bool debug,
                       ROI_UPDATE_MODE roi_update_mode = ROI_UPDATE_MODE::SKIP_ALL,
                       ACCUMULATION_MODE accumulation_mode = ACCUMULATION_MODE::ON_THE_FLY,
                       bool save_disps_steps = false) : 
                                     imgs(imgs),
                                     roi(roi),
                                     scalefactor(scalefactor),
                                     interp_type(interp_type),
                                     subregion_type(subregion_type), 
                                     r(r), 
                                     num_threads(num_threads), 
                                     cutoff_corrcoef(cutoff_corrcoef), 
                                     update_corrcoef(update_corrcoef),
                                     prctile_corrcoef(prctile_corrcoef), 
                                     roi_update_mode(roi_update_mode),
                                     accumulation_mode(accumulation_mode),
                                     save_disps_steps(save_disps_steps || debug),  // Enforce if debug is true
                                     debug(debug) { }    
    
    // Constructor with config - can optionally override roi_update_mode and accumulation_mode
    DIC_analysis_input(const std::vector<Image2D>&, const ROI2D&, difference_type, INTERP, SUBREGION, difference_type, difference_type, DIC_analysis_config, bool, 
                       ROI_UPDATE_MODE roi_update_mode_override = static_cast<ROI_UPDATE_MODE>(-1),
                       ACCUMULATION_MODE accumulation_mode_override = static_cast<ACCUMULATION_MODE>(-1),
                       bool save_disps_steps_override = false);
        
    // Static factory methods ------------------------------------------------//
    static DIC_analysis_input load(std::ifstream&);
    static DIC_analysis_input load(const std::string&);
            
    // Interface functions ---------------------------------------------------//
    friend void save(const DIC_analysis_input&, std::ofstream&);
    friend void save(const DIC_analysis_input&, const std::string&);
        
    std::vector<Image2D> imgs;
    ROI2D roi;
    difference_type scalefactor;
    INTERP interp_type;
    SUBREGION subregion_type; 
    difference_type r; 
    difference_type num_threads;
    double cutoff_corrcoef;
    double update_corrcoef;
    double prctile_corrcoef;
    ROI_UPDATE_MODE roi_update_mode;  // Mode for ROI boundary update (SKIP_ALL or SKIP_INVALID)
    ACCUMULATION_MODE accumulation_mode;  // Mode for displacement accumulation (ON_THE_FLY or POST_PROCESS)
    bool save_disps_steps;  // If true, save step displacements data for debugging (enforced when debug=true)
    bool debug;
};
 
enum class PERSPECTIVE { EULERIAN, LAGRANGIAN };
 
struct DIC_analysis_output final {
    typedef ROI2D::difference_type                              difference_type;
        
    // Rule of 5 and destructor ----------------------------------------------//    
    DIC_analysis_output() : units_per_pixel() { }
    DIC_analysis_output(const DIC_analysis_output&) = default;
    DIC_analysis_output(DIC_analysis_output&&) = default;
    DIC_analysis_output& operator=(const DIC_analysis_output&) = default;
    DIC_analysis_output& operator=(DIC_analysis_output&&) = default;
    ~DIC_analysis_output() noexcept = default;
    
    // Additional constructors -----------------------------------------------//    
    DIC_analysis_output(const std::vector<Disp2D> &disps, PERSPECTIVE perspective_type, const std::string &units, double units_per_pixel) : 
        disps(disps), perspective_type(perspective_type), units(units), units_per_pixel(units_per_pixel) { }
    
    // Static factory methods ------------------------------------------------//
    static DIC_analysis_output load(std::ifstream&);
    static DIC_analysis_output load(const std::string&);
            
    // Interface functions ---------------------------------------------------//
    friend void save(const DIC_analysis_output&, std::ofstream&);
    friend void save(const DIC_analysis_output&, const std::string&);
        
    std::vector<Disp2D> disps;
    PERSPECTIVE perspective_type;
    std::string units;
    double units_per_pixel;
};
 
// Step displacement data for debugging and analysis
struct DIC_analysis_step_data final {
    typedef ROI2D::difference_type                              difference_type;
    
    // Rule of 5 and destructor
    DIC_analysis_step_data() = default;
    DIC_analysis_step_data(const DIC_analysis_step_data&) = default;
    DIC_analysis_step_data(DIC_analysis_step_data&&) = default;
    DIC_analysis_step_data& operator=(const DIC_analysis_step_data&) = default;
    DIC_analysis_step_data& operator=(DIC_analysis_step_data&&) = default;
    ~DIC_analysis_step_data() noexcept = default;
    
    // Static factory methods
    static DIC_analysis_step_data load(std::ifstream&);
    static DIC_analysis_step_data load(const std::string&);
    
    // Interface functions
    friend void save(const DIC_analysis_step_data&, std::ofstream&);
    friend void save(const DIC_analysis_step_data&, const std::string&);
    
    // Step displacements: step_disps[i] = displacement from frame step_ref_idx[i] to frame i+1
    std::vector<Disp2D> step_disps;
    // ROIs used when computing each step displacement
    std::vector<ROI2D> step_rois;
    // Reference frame index for each step
    std::vector<difference_type> step_ref_idx;
};
 
DIC_analysis_output DIC_analysis(const DIC_analysis_input&);
DIC_analysis_output DIC_analysis_sequential(const DIC_analysis_input&);
 
// Seed-based parallel DIC analysis ------------------------------------------//
struct SeedParams final {
    typedef ROI2D::difference_type                              difference_type;
    
    // Seed position in the reference image
    difference_type x;
    difference_type y;
    
    // Displacement and deformation gradient (9 parameters total)
    double u;          // displacement in x
    double v;          // displacement in y
    double du_dx;      // deformation gradient
    double du_dy;
    double dv_dx;
    double dv_dy;
    double corrcoef;   // correlation coefficient (quality metric)
    
    // Constructors
    SeedParams() : x(0), y(0), u(0.0), v(0.0), du_dx(0.0), du_dy(0.0), dv_dx(0.0), dv_dy(0.0), corrcoef(0.0) { }
    SeedParams(difference_type x, difference_type y) : 
        x(x), y(y), u(0.0), v(0.0), du_dx(0.0), du_dy(0.0), dv_dx(0.0), dv_dy(0.0), corrcoef(0.0) { }
    
    // Convert to Array2D format compatible with RGDIC functions
    // params = {p1, p2, v, u, dv_dp1, dv_dp2, du_dp1, du_dp2, corr_coef, diff_norm}
    // Note: p1=y (row), p2=x (col), and gradients map: dv_dp1=dv_dy, dv_dp2=dv_dx, du_dp1=du_dy, du_dp2=du_dx
    Array2D<double> to_array(double diff_norm = 0.0) const {
        Array2D<double> params(10, 1);
        params(0) = static_cast<double>(y);      // p1 (row)
        params(1) = static_cast<double>(x);      // p2 (col)
        params(2) = v;                            // v displacement
        params(3) = u;                            // u displacement
        params(4) = dv_dy;                        // dv_dp1
        params(5) = dv_dx;                        // dv_dp2
        params(6) = du_dy;                        // du_dp1
        params(7) = du_dx;                        // du_dp2
        params(8) = corrcoef;                     // correlation coefficient
        params(9) = diff_norm;                    // difference norm (convergence metric)
        return params;
    }
    
    // Create SeedParams from Array2D format
    static SeedParams from_array(const Array2D<double>& params) {
        SeedParams seed;
        seed.y = static_cast<difference_type>(params(0));  // p1 (row)
        seed.x = static_cast<difference_type>(params(1));  // p2 (col)
        seed.v = params(2);
        seed.u = params(3);
        seed.dv_dy = params(4);  // dv_dp1
        seed.dv_dx = params(5);  // dv_dp2
        seed.du_dy = params(6);  // du_dp1
        seed.du_dx = params(7);  // du_dp2
        seed.corrcoef = params(8);
        return seed;
    }
};
 
struct SeedConvergence final {
    int num_iterations;
    double diffnorm;
    
    SeedConvergence() : num_iterations(0), diffnorm(0.0) { }
    SeedConvergence(int iter, double norm) : num_iterations(iter), diffnorm(norm) { }
};
 
struct SeedAnalysisResult final {
    std::vector<SeedParams> seeds;
    std::vector<SeedConvergence> convergence;
    bool success;
    
    SeedAnalysisResult() : success(false) { }
};
 
// Lightweight seed analysis for predicting reference updates
SeedAnalysisResult analyze_seeds(
    const details::subregion_nloptimizer &sr_nloptimizer,
    const Array2D<double>& ref_gs,
    const ROI2D& roi,
    const std::vector<SeedParams>& seed_positions,
    ROI2D::difference_type radius,
    int cutoff_iteration,
    double cutoff_max_diffnorm,
    double cutoff_max_corrcoef,
    bool debug = true
);
 
// Structure to hold precomputed seed data for each frame
struct SeedComputationData final {
    ROI2D roi;                    // Updated ROI for this frame
    std::vector<SeedParams> seed_params_by_region;  // One seed param per region
    details::subregion_nloptimizer sr_nloptimizer;
    
    SeedComputationData() = default;
    SeedComputationData(const ROI2D& roi, const std::vector<SeedParams>& params, const details::subregion_nloptimizer &sr_nloptimizer) : 
        roi(roi), seed_params_by_region(params), sr_nloptimizer(sr_nloptimizer) { }
};
 
// Compute seed parameters for all frames with ROI updates
// Returns precomputed data (roi, seed_params) for each successfully analyzed frame
std::vector<SeedComputationData> compute_only_seed_points(
    const Array2D<double>& A_ref,
    const std::vector<Array2D<double>>& A_curs,
    const ROI2D& roi,
    ROI2D::difference_type scalefactor,
    INTERP interp_type,
    SUBREGION subregion_type,
    const std::vector<SeedParams>& seeds_by_region,  // One seed per region
    double cutoff_corrcoef,
    ROI2D::difference_type region_idx = 0,
    bool debug = false
);
 
// Compute displacements using precomputed seed parameters
Disp2D compute_displacements(
    const details::subregion_nloptimizer &sr_nloptimizer,
    const ROI2D& roi_reduced,
    const SeedParams& seedparams,
    ROI2D::difference_type scalefactor,
    double cutoff_corrcoef,
    ROI2D::difference_type region_idx,
    bool debug
) ;
 
// Update seed positions based on displacement (seed propagation)
std::vector<SeedParams> propagate_seeds(
    const std::vector<SeedParams>& seeds,
    ROI2D::difference_type spacing
);
 
struct DIC_analysis_parallel_input final {
    typedef ROI2D::difference_type                              difference_type;
    
    // Base DIC input parameters
    DIC_analysis_input base_input;
    
    // Seed parameters (one per region - provided by user, not generated)
    std::vector<SeedParams> seeds_by_region;
    
    // Whether seeds are already optimized (skip optimization step if true)
    bool seeds_are_optimized;
    
    // Seed-based parallelization parameters
    double cutoff_max_diffnorm;          // Diffnorm threshold for failure prediction
    double cutoff_max_corrcoef;          // Corrcoef threshold for failure prediction
    
    // Constructors
    DIC_analysis_parallel_input() : seeds_are_optimized(false), cutoff_max_diffnorm(0.1), cutoff_max_corrcoef(0.5) { }
    
    DIC_analysis_parallel_input(const DIC_analysis_input& base, const std::vector<SeedParams>& seeds, bool optimized = false) : 
        base_input(base), seeds_by_region(seeds), seeds_are_optimized(optimized),
        cutoff_max_diffnorm(0.1), cutoff_max_corrcoef(0.5) { }
};
 
// Parallel DIC analysis using seed-based failure prediction
DIC_analysis_output DIC_analysis_parallel(const DIC_analysis_parallel_input&);
 
// Matlab-ABR-style DIC analysis using a fixed reference, manual/preset seeds,
// seed-analysis first, then DIC evaluation.
DIC_analysis_output matlab_DIC_analysis_parallel(const DIC_analysis_parallel_input&);
 
// Sequential DIC analysis with user-provided seeds
DIC_analysis_output DIC_analysis_sequential(
    const DIC_analysis_input& DIC_input,
    const std::vector<SeedParams>& seeds_by_region = {},
    bool seeds_are_optimized = false
);
 
DIC_analysis_output DIC_analysis_sequential(const DIC_analysis_parallel_input&);
 
// Sequential Matlab-ABR-style DIC analysis with user-provided seeds.
DIC_analysis_output matlab_DIC_analysis_sequential(
    const DIC_analysis_input& DIC_input,
    const std::vector<SeedParams>& seeds_by_region = {},
    bool seeds_are_optimized = false
);
 
DIC_analysis_output matlab_DIC_analysis_sequential(const DIC_analysis_parallel_input&);
 
// Exact MATLAB-ncorr-mirroring DIC analysis. Same seed/segment pipeline as
// matlab_DIC_analysis_*, but the multi-reference chain composition uses
// exact_add_with_rois(), which mirrors ncorr_alg_addanalysis + ncorr_alg_adddisp:
//  - per-region extrapolation via expand_filt (4-neighbor mean iteration)
//  - biquintic B-spline interpolation on the extrapolated plot
//  - reduced-grid walk with strict in_nlinfo mask check
// This fixes the chain-induced jump at the first segment boundary that
// add_with_rois() exhibited (see seam diagnostics in matlab_DIC_analysis_*).
DIC_analysis_output exact_matlab_DIC_analysis_parallel(const DIC_analysis_parallel_input&);
DIC_analysis_output exact_matlab_DIC_analysis_sequential(
    const DIC_analysis_input& DIC_input,
    const std::vector<SeedParams>& seeds_by_region = {},
    bool seeds_are_optimized = false
);
DIC_analysis_output exact_matlab_DIC_analysis_sequential(const DIC_analysis_parallel_input&);
 
// MATLAB-style chain composition: extrapdata + biquintic + reduced-grid walk.
Disp2D exact_add_with_rois(const std::vector<Disp2D>& disps,
                           const std::vector<ROI2D>& rois);
 
// RGDIC with user-provided seeds (single frame pair)
Disp2D RGDIC_with_seeds(const Array2D<double>& A_ref, 
                        const Array2D<double>& A_cur, 
                        const ROI2D& roi,
                        const DIC_analysis_parallel_input& input);
 
// Conversion between Lagrangian and Eulerian displacements ------------------//
DIC_analysis_output change_perspective(const DIC_analysis_output&, INTERP);
 
// Conversion with Matlab-compatible sign inversion for Eulerian perspective -//
DIC_analysis_output change_perspective_with_inversion(const DIC_analysis_output&, INTERP);
 
// set units -----------------------------------------------------------------//
DIC_analysis_output set_units(const DIC_analysis_output&, const std::string&, double);
 
// Correlation filtering -----------------------------------------------------//
DIC_analysis_output filter_by_correlation(const DIC_analysis_output&, double);
 
// strain_analysis -----------------------------------------------------------//
Strain2D LS_strain(const Disp2D&, PERSPECTIVE, double, SUBREGION, ROI2D::difference_type); 
 
struct strain_analysis_input final {
    typedef ROI2D::difference_type                              difference_type;
        
    // Rule of 5 and destructor ----------------------------------------------//    
    strain_analysis_input() : r() { }
    strain_analysis_input(const strain_analysis_input&) = default;
    strain_analysis_input(strain_analysis_input&&) = default;
    strain_analysis_input& operator=(const strain_analysis_input&) = default;
    strain_analysis_input& operator=(strain_analysis_input&&) = default;
    ~strain_analysis_input() noexcept = default;
    
    // Additional constructors -----------------------------------------------//    
    strain_analysis_input(const DIC_analysis_input &DIC_input,
                          const DIC_analysis_output &DIC_output,
                          SUBREGION subregion_type,
                          difference_type r) : DIC_input(DIC_input),
                                               DIC_output(DIC_output),
                                               subregion_type(subregion_type),
                                               r(r){ }   
    
    // Static factory methods ------------------------------------------------//
    static strain_analysis_input load(std::ifstream&);
    static strain_analysis_input load(const std::string&);
            
    // Interface functions ---------------------------------------------------//
    friend void save(const strain_analysis_input&, std::ofstream&);
    friend void save(const strain_analysis_input&, const std::string&);
        
    DIC_analysis_input DIC_input;
    DIC_analysis_output DIC_output;
    SUBREGION subregion_type; 
    difference_type r;
};
 
struct strain_analysis_output final {
    typedef ROI2D::difference_type                              difference_type;
        
    // Rule of 5 and destructor ----------------------------------------------//    
    strain_analysis_output() = default;
    strain_analysis_output(const strain_analysis_output&) = default;
    strain_analysis_output(strain_analysis_output&&) = default;
    strain_analysis_output& operator=(const strain_analysis_output&) = default;
    strain_analysis_output& operator=(strain_analysis_output&&) = default;
    ~strain_analysis_output() noexcept = default;
    
    // Additional constructors -----------------------------------------------//    
    strain_analysis_output(const std::vector<Strain2D> &strains) : strains(strains) { }
    
    // Static factory methods ------------------------------------------------//
    static strain_analysis_output load(std::ifstream&);
    static strain_analysis_output load(const std::string&);
            
    // Interface functions ---------------------------------------------------//
    friend void save(const strain_analysis_output&, std::ofstream&);
    friend void save(const strain_analysis_output&, const std::string&);
        
    std::vector<Strain2D> strains;
};
 
strain_analysis_output strain_analysis(const strain_analysis_input&);
 
// Interface functions for viewing and saving ncorr related data -------------//
void imshow_ncorr_data_over_img(const Image2D&, const Data2D&, ROI2D::difference_type = -1);
 
void save_ncorr_data_over_img(const std::string&, 
                              const Image2D&, 
                              const Data2D&, 
                              double, 
                              double, 
                              double, 
                              bool, 
                              bool, 
                              bool, 
                              const std::string&, 
                              double, 
                              double, 
                              double,
                              ROI2D::difference_type,
                              int);
 
void save_ncorr_data_over_img_video(const std::string&, 
                                    const std::vector<Image2D>&, 
                                    const std::vector<Data2D>&, 
                                    double, 
                                    double, 
                                    double, 
                                    double, 
                                    bool, 
                                    bool, 
                                    bool, 
                                    const std::string&, 
                                    double, 
                                    double, 
                                    double,
                                    ROI2D::difference_type,
                                    int, 
                                    double, 
                                    int);
 
enum class DISP { U, V };
void save_DIC_video(const std::string&, 
                    const DIC_analysis_input&, 
                    const DIC_analysis_output&, 
                    DISP, 
                    double, 
                    double, 
                    double = std::numeric_limits<double>::quiet_NaN(), 
                    double = std::numeric_limits<double>::quiet_NaN(), 
                    bool = true, 
                    bool = true, 
                    bool = true, 
                    double = -1.0,
                    double = 1.0, 
                    ROI2D::difference_type = 11,
                    int = cv::COLORMAP_JET,
                    double = 2.0, 
                    int = cv::VideoWriter::fourcc('M','J','P','G'));
 
enum class STRAIN { EYY, EXY, EXX };
void save_strain_video(const std::string&, 
                       const strain_analysis_input&, 
                       const strain_analysis_output&, 
                       STRAIN, 
                       double, 
                       double, 
                       double = std::numeric_limits<double>::quiet_NaN(), 
                       double = std::numeric_limits<double>::quiet_NaN(), 
                       bool = true, 
                       bool = true, 
                       bool = true, 
                       double = -1.0,
                       double = 1.0, 
                       ROI2D::difference_type = 11,
                       int = cv::COLORMAP_JET,
                       double = 2.0, 
                       int = cv::VideoWriter::fourcc('M','J','P','G'));
// ---------------------------------------------------------------------------//
 
}
 
#endif  /* NCORR_H */