CPPxDIC
C++ Digital Image Correlation library for fingertip 3D reconstruction
Loading...
Searching...
No Matches
dic_analysis.cpp
Go to the documentation of this file.
1
5#include "dic_analysis.h"
6#include "mat_reader.h"
7#include "utils.h"
8#include "step_d_workflow.h"
9#include "data_serializer.h"
10#include "strain_computation.h"
11#include "surface_stitching.h"
12#include "temporal_filter.h"
13#include "face_isotropy.h"
14#include "visualization.h"
15#include "logging.h"
16#include <iostream>
17#include <chrono>
18#include <filesystem>
19#include <sstream>
20#include <iomanip>
21#include <matio.h>
22#include <limits>
23#include <Eigen/Dense>
24#include <cctype>
25#include <cmath>
26#include <algorithm>
27#include <array>
28#include <numeric>
29#include <set>
30
31using namespace ncorr;
32using namespace cppxdic;
33
34DicAnalysis::DicAnalysis(const Config& config) : config_(config) {
35 // Library entry point: make sure the configured stereopair table is active even
36 // when the caller is not main.cpp (tests, tools, other drivers).
38}
39
40bool DicAnalysis::dicDeformationAnalysis(const std::vector<int>& trial_target) {
41 LOG_INFO << "Starting Deformation/Strain Analysis (Step F)...";
42 LOG_INFO << "NOTE: This requires DIC3Dcombined from Step E";
43
44 try {
45 for (int trial : trial_target) {
46 // Step F works on combined 3D reconstruction, not per-pair
47 std::string output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
48
49 // Load DIC3Dcombined from Step E output using configured format
50 auto serializer = cppxdic::DataSerializer::create(config_.data_format);
51 std::string dic3d_file = Utils::buildDic3DCombinedFilePath(output_dir, config_.num_pair, serializer->extension());
52
53 if (!std::filesystem::exists(dic3d_file)) {
54 LOG_ERROR << "DIC3Dcombined file not found: " << dic3d_file;
55 LOG_ERROR << "You must run Step E (dic3DReconstruction) first!";
56 return false;
57 }
58
59 LOG_INFO << "Loading DIC3Dcombined from: " << dic3d_file;
60
61 DIC3Dcombined dic3d;
62 if (!serializer->loadDIC3Dcombined(dic3d_file, dic3d)) {
63 LOG_ERROR << "Failed to load DIC3Dcombined structure";
64 return false;
65 }
66
67 if (dic3d.Points3D.empty() || dic3d.Faces.empty()) {
68 LOG_ERROR << "DIC3Dcombined has no 3D data";
69 return false;
70 }
71
72 LOG_INFO << " Loaded: " << dic3d.Points3D.size() << " frames, "
73 << dic3d.Points3D[0].x.size() << " points, "
74 << dic3d.Faces.size() / 3 << " faces";
75 LOG_INFO << " corrComb: " << dic3d.corrComb.size() << " frames"
76 << (dic3d.corrComb.empty() ? "" : " (" + std::to_string(dic3d.corrComb[0].size()) + " per frame)");
77 LOG_INFO << " FaceCorrComb: " << dic3d.FaceCorrComb.size() << " frames";
78 LOG_INFO << " FaceCentroids: " << dic3d.FaceCentroids.size() << " frames";
79 LOG_INFO << " Disp.DispVec: " << dic3d.Disp.DispVec.size() << " frames";
80 LOG_INFO << " Disp.DispMgn: " << dic3d.Disp.DispMgn.size() << " frames";
81
82 // Diagnostic: check for NaN/Inf in Points3D and face validity
83 {
84 size_t inf_pts_total = 0;
85 size_t bad_face_total = 0;
86 size_t nP = dic3d.Points3D.empty() ? 0 : dic3d.Points3D[0].x.size();
87 for (size_t f = 0; f < std::min(dic3d.Points3D.size(), size_t(3)); ++f) {
88 size_t nan_pts = 0;
89 for (size_t i = 0; i < nP; ++i) {
90 if (std::isnan(dic3d.Points3D[f].x[i]) || std::isnan(dic3d.Points3D[f].y[i]) || std::isnan(dic3d.Points3D[f].z[i]))
91 nan_pts++;
92 if (std::isinf(dic3d.Points3D[f].x[i]) || std::isinf(dic3d.Points3D[f].y[i]) || std::isinf(dic3d.Points3D[f].z[i]))
93 inf_pts_total++;
94 }
95 LOG_DEBUG << " [diag] Frame " << f << ": " << nan_pts << "/" << nP << " NaN points";
96 }
97 size_t nF = dic3d.Faces.size() / 3;
98 for (size_t i = 0; i < nF; ++i) {
99 int v0 = dic3d.Faces[i*3], v1 = dic3d.Faces[i*3+1], v2 = dic3d.Faces[i*3+2];
100 if (v0 < 0 || v1 < 0 || v2 < 0 || (size_t)v0 >= nP || (size_t)v1 >= nP || (size_t)v2 >= nP)
101 bad_face_total++;
102 }
103 if (bad_face_total > 0)
104 LOG_WARN << " [diag] " << bad_face_total << "/" << nF << " faces have out-of-bounds vertex indices!";
105 if (inf_pts_total > 0)
106 LOG_WARN << " [diag] " << inf_pts_total << " Inf values in Points3D!";
107 }
108
109 // Convert Points3D to Eigen::Vector3d format for deformation computation
110 LOG_INFO << "\nConverting data to Eigen format...";
111
112 // Reference frame (frame 0)
113 std::vector<Eigen::Vector3d> vertices_ref;
114 size_t nPoints = dic3d.Points3D[0].x.size();
115 size_t nFaces = dic3d.Faces.size() / 3;
116 vertices_ref.reserve(nPoints);
117 for (size_t i = 0; i < nPoints; ++i) {
118 vertices_ref.emplace_back(
119 dic3d.Points3D[0].x[i],
120 dic3d.Points3D[0].y[i],
121 dic3d.Points3D[0].z[i]
122 );
123 }
124
125 // All frames
126 std::vector<std::vector<Eigen::Vector3d>> vertices_all_frames;
127 vertices_all_frames.reserve(dic3d.Points3D.size());
128 for (const auto& frame_pts : dic3d.Points3D) {
129 std::vector<Eigen::Vector3d> frame_verts;
130 frame_verts.reserve(frame_pts.x.size());
131 for (size_t i = 0; i < frame_pts.x.size(); ++i) {
132 frame_verts.emplace_back(
133 frame_pts.x[i],
134 frame_pts.y[i],
135 frame_pts.z[i]
136 );
137 }
138 vertices_all_frames.push_back(std::move(frame_verts));
139 }
140
141 LOG_INFO << " Converted " << vertices_all_frames.size() << " frames";
142
143 auto count_nan_points = [](const std::vector<std::vector<Eigen::Vector3d>>& frames) {
144 std::vector<size_t> counts;
145 counts.reserve(frames.size());
146 for (const auto& frame : frames) {
147 size_t nan_count = 0;
148 for (const auto& p : frame) {
149 if (p.hasNaN()) {
150 nan_count++;
151 }
152 }
153 counts.push_back(nan_count);
154 }
155 return counts;
156 };
157
158 auto count_points_with_any_nan = [](const std::vector<std::vector<Eigen::Vector3d>>& frames) {
159 if (frames.empty()) {
160 return size_t(0);
161 }
162 size_t n_points_local = frames.front().size();
163 size_t bad_points = 0;
164 for (size_t ipt = 0; ipt < n_points_local; ++ipt) {
165 bool any_nan = false;
166 for (const auto& frame : frames) {
167 if (ipt >= frame.size() || frame[ipt].hasNaN()) {
168 any_nan = true;
169 break;
170 }
171 }
172 if (any_nan) {
173 bad_points++;
174 }
175 }
176 return bad_points;
177 };
178
179 auto count_valid_faces = [&](const std::vector<Eigen::Vector3d>& verts) {
180 size_t valid_faces = 0;
181 size_t total_faces_local = dic3d.Faces.size() / 3;
182 for (size_t iface = 0; iface < total_faces_local; ++iface) {
183 int v0 = dic3d.Faces[iface * 3 + 0];
184 int v1 = dic3d.Faces[iface * 3 + 1];
185 int v2 = dic3d.Faces[iface * 3 + 2];
186 if (v0 < 0 || v1 < 0 || v2 < 0) {
187 continue;
188 }
189 if (static_cast<size_t>(v0) >= verts.size() ||
190 static_cast<size_t>(v1) >= verts.size() ||
191 static_cast<size_t>(v2) >= verts.size()) {
192 continue;
193 }
194 if (verts[v0].hasNaN() || verts[v1].hasNaN() || verts[v2].hasNaN()) {
195 continue;
196 }
197 valid_faces++;
198 }
199 return valid_faces;
200 };
201
202 auto compute_max_displacement = [](const std::vector<std::vector<Eigen::Vector3d>>& frames,
203 const std::vector<Eigen::Vector3d>& ref) {
204 double max_disp = 0.0;
205 for (const auto& frame : frames) {
206 size_t n_local = std::min(frame.size(), ref.size());
207 for (size_t ipt = 0; ipt < n_local; ++ipt) {
208 if (!frame[ipt].allFinite() || !ref[ipt].allFinite()) {
209 continue;
210 }
211 max_disp = std::max(max_disp, (frame[ipt] - ref[ipt]).norm());
212 }
213 }
214 return max_disp;
215 };
216
217 const size_t raw_any_nan_points = count_points_with_any_nan(vertices_all_frames);
218 const auto raw_nan_counts = count_nan_points(vertices_all_frames);
219 const double raw_max_disp = compute_max_displacement(vertices_all_frames, vertices_ref);
220 LOG_DEBUG << " Raw points with any NaN over time: " << raw_any_nan_points
221 << "/" << nPoints;
222 if (!raw_nan_counts.empty()) {
223 LOG_DEBUG << " Raw NaN points in frame 0: " << raw_nan_counts[0]
224 << "/" << nPoints;
225 }
226 LOG_DEBUG << " Raw valid faces in frame 0: "
227 << count_valid_faces(vertices_all_frames[0])
228 << "/" << nFaces;
229 LOG_DEBUG << " Raw max displacement magnitude: " << raw_max_disp;
230
231 // Apply temporal filtering to displacement fields
232 bool used_temporal_filtering = false;
233 if (config_.step_f_temporal_filtering && vertices_all_frames.size() > 3) {
234 LOG_INFO << "\nApplying temporal filtering...";
235
236 // Organize data for filtering: nPoints x nFrames
237 size_t nFrames = vertices_all_frames.size();
238 std::vector<std::vector<double>> disp_x(nPoints, std::vector<double>(nFrames));
239 std::vector<std::vector<double>> disp_y(nPoints, std::vector<double>(nFrames));
240 std::vector<std::vector<double>> disp_z(nPoints, std::vector<double>(nFrames));
241
242 for (size_t iframe = 0; iframe < nFrames; ++iframe) {
243 for (size_t ipt = 0; ipt < nPoints; ++ipt) {
244 disp_x[ipt][iframe] = vertices_all_frames[iframe][ipt].x() - vertices_ref[ipt].x();
245 disp_y[ipt][iframe] = vertices_all_frames[iframe][ipt].y() - vertices_ref[ipt].y();
246 disp_z[ipt][iframe] = vertices_all_frames[iframe][ipt].z() - vertices_ref[ipt].z();
247 }
248 }
249
250 // Filter displacement components
251 double freq_filt = config_.step_f_freq_filt;
252 double freq_acq = config_.vid_sample_freq;
253 auto [filt_x, filt_y, filt_z] = filterTime3D(disp_x, disp_y, disp_z, freq_filt, freq_acq);
254
255 // Reconstruct filtered vertex positions
256 auto vertices_filtered = vertices_all_frames;
257 for (size_t iframe = 0; iframe < nFrames; ++iframe) {
258 for (size_t ipt = 0; ipt < nPoints; ++ipt) {
259 vertices_filtered[iframe][ipt].x() = vertices_ref[ipt].x() + filt_x[ipt][iframe];
260 vertices_filtered[iframe][ipt].y() = vertices_ref[ipt].y() + filt_y[ipt][iframe];
261 vertices_filtered[iframe][ipt].z() = vertices_ref[ipt].z() + filt_z[ipt][iframe];
262 }
263 }
264
265 const size_t filtered_any_nan_points = count_points_with_any_nan(vertices_filtered);
266 const auto filtered_nan_counts = count_nan_points(vertices_filtered);
267 const size_t filtered_valid_faces = count_valid_faces(vertices_filtered[0]);
268 const double filtered_max_disp = compute_max_displacement(vertices_filtered, vertices_ref);
269 const bool filtered_exploded =
270 !std::isfinite(filtered_max_disp) ||
271 filtered_max_disp > std::max(1.0, raw_max_disp) * 100.0;
272
273 LOG_DEBUG << " Filtered points with any NaN over time: " << filtered_any_nan_points
274 << "/" << nPoints;
275 if (!filtered_nan_counts.empty()) {
276 LOG_DEBUG << " Filtered NaN points in frame 0: " << filtered_nan_counts[0]
277 << "/" << nPoints;
278 }
279 LOG_DEBUG << " Filtered valid faces in frame 0: "
280 << filtered_valid_faces << "/" << nFaces;
281 LOG_DEBUG << " Filtered max displacement magnitude: " << filtered_max_disp;
282
283 // MATLAB's filter only keeps point tracks that are valid for the whole time series.
284 // On stitched C++ reconstructions this can erase nearly the entire mesh, leaving
285 // TCPE with no valid triangles. The current C++ filter can also become numerically
286 // unstable and amplify otherwise small motions into absurd coordinates. In either
287 // case we keep the raw geometry instead of feeding corrupted data into Step F.
288 if (filtered_valid_faces == 0 || filtered_valid_faces * 20 < nFaces || filtered_exploded) {
289 LOG_WARN << " Temporal filtering produced unusable geometry; "
290 << "keeping unfiltered geometry for Step F";
291 } else {
292 vertices_all_frames = std::move(vertices_filtered);
293 for (size_t ipt = 0; ipt < nPoints; ++ipt) {
294 vertices_ref[ipt] = vertices_all_frames[0][ipt];
295 }
296 used_temporal_filtering = true;
297 LOG_INFO << " ✓ Temporal filtering applied (freqFilt=" << freq_filt
298 << " Hz, freqAcq=" << freq_acq << " Hz)";
299 }
300 } else if (config_.step_f_temporal_filtering) {
301 LOG_INFO << "\nSkipping temporal filtering (too few frames: "
302 << vertices_all_frames.size() << ")";
303 } else {
304 LOG_INFO << "\nTemporal filtering disabled (step_f_temporal_filtering=false)";
305 }
306
307 // Update Points3D with filtered data (matching MATLAB line 115)
308 LOG_INFO << "\nUpdating Points3D with filtered data...";
309 dic3d.Points3D.clear();
310 dic3d.Points3D.resize(vertices_all_frames.size());
311 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
312 Point3DFrame& frame_pts = dic3d.Points3D[iframe];
313 frame_pts.x.resize(nPoints);
314 frame_pts.y.resize(nPoints);
315 frame_pts.z.resize(nPoints);
316 for (size_t ipt = 0; ipt < nPoints; ++ipt) {
317 frame_pts.x[ipt] = vertices_all_frames[iframe][ipt].x();
318 frame_pts.y[ipt] = vertices_all_frames[iframe][ipt].y();
319 frame_pts.z[ipt] = vertices_all_frames[iframe][ipt].z();
320 }
321 }
322 if (used_temporal_filtering) {
323 LOG_INFO << " ✓ Points3D updated with filtered data";
324 } else {
325 LOG_INFO << " ✓ Points3D kept from unfiltered geometry";
326 }
327
328 // Recompute displacement based on filtered Points3D (matching MATLAB lines 69-71)
329 LOG_INFO << "\nRecomputing displacement after filtering...";
330 dic3d.Disp.DispVec.clear();
331 dic3d.Disp.DispMgn.clear();
332 dic3d.Disp.DispVec.resize(vertices_all_frames.size());
333 dic3d.Disp.DispMgn.resize(vertices_all_frames.size());
334
335 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
336 std::vector<double>& disp_vec = dic3d.Disp.DispVec[iframe];
337 std::vector<double>& disp_mgn = dic3d.Disp.DispMgn[iframe];
338
339 disp_vec.resize(nPoints * 3);
340 disp_mgn.resize(nPoints);
341
342 for (size_t ipt = 0; ipt < nPoints; ++ipt) {
343 double dx = vertices_all_frames[iframe][ipt].x() - vertices_ref[ipt].x();
344 double dy = vertices_all_frames[iframe][ipt].y() - vertices_ref[ipt].y();
345 double dz = vertices_all_frames[iframe][ipt].z() - vertices_ref[ipt].z();
346
347 disp_vec[ipt * 3 + 0] = dx;
348 disp_vec[ipt * 3 + 1] = dy;
349 disp_vec[ipt * 3 + 2] = dz;
350
351 disp_mgn[ipt] = std::sqrt(dx*dx + dy*dy + dz*dz);
352 }
353 }
354 LOG_INFO << " ✓ Displacement recomputed for " << vertices_all_frames.size() << " frames";
355
356 // Recompute face centroids based on filtered Points3D (matching MATLAB lines 64-66)
357 LOG_INFO << "\nRecomputing face centroids after filtering...";
358 dic3d.FaceCentroids.clear();
359 dic3d.FaceCentroids.resize(vertices_all_frames.size());
360
361 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
362 std::vector<double>& centroids = dic3d.FaceCentroids[iframe];
363 centroids.resize(nFaces * 3);
364
365 size_t nVerts = vertices_all_frames[iframe].size();
366 for (size_t iface = 0; iface < nFaces; ++iface) {
367 int v0 = dic3d.Faces[iface * 3 + 0];
368 int v1 = dic3d.Faces[iface * 3 + 1];
369 int v2 = dic3d.Faces[iface * 3 + 2];
370
371 if (v0 >= 0 && static_cast<size_t>(v0) < nVerts &&
372 v1 >= 0 && static_cast<size_t>(v1) < nVerts &&
373 v2 >= 0 && static_cast<size_t>(v2) < nVerts) {
374 double cx = (vertices_all_frames[iframe][v0].x() +
375 vertices_all_frames[iframe][v1].x() +
376 vertices_all_frames[iframe][v2].x()) / 3.0;
377 double cy = (vertices_all_frames[iframe][v0].y() +
378 vertices_all_frames[iframe][v1].y() +
379 vertices_all_frames[iframe][v2].y()) / 3.0;
380 double cz = (vertices_all_frames[iframe][v0].z() +
381 vertices_all_frames[iframe][v1].z() +
382 vertices_all_frames[iframe][v2].z()) / 3.0;
383
384 centroids[iface * 3 + 0] = cx;
385 centroids[iface * 3 + 1] = cy;
386 centroids[iface * 3 + 2] = cz;
387 } else {
388 centroids[iface * 3 + 0] = std::numeric_limits<double>::quiet_NaN();
389 centroids[iface * 3 + 1] = std::numeric_limits<double>::quiet_NaN();
390 centroids[iface * 3 + 2] = std::numeric_limits<double>::quiet_NaN();
391 }
392 }
393 }
394 LOG_INFO << " ✓ Face centroids recomputed for " << vertices_all_frames.size() << " frames";
395
396 // Recompute face correlation (worst of 3 vertices) (matching MATLAB line 61)
397 // MATLAB: DIC3D.FaceCorrComb{ii} = max(DIC3D.corrComb{ii}(F), [], 2);
398 LOG_INFO << "\nRecomputing face correlation after filtering...";
399 if (!dic3d.corrComb.empty() && dic3d.corrComb.size() >= vertices_all_frames.size()) {
400 dic3d.FaceCorrComb.clear();
401 dic3d.FaceCorrComb.resize(vertices_all_frames.size());
402
403 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
404 std::vector<double>& face_corr = dic3d.FaceCorrComb[iframe];
405 face_corr.resize(nFaces);
406
407 const auto& point_corr = dic3d.corrComb[iframe];
408 size_t nCorr = point_corr.size();
409
410 for (size_t iface = 0; iface < nFaces; ++iface) {
411 int v0 = dic3d.Faces[iface * 3 + 0];
412 int v1 = dic3d.Faces[iface * 3 + 1];
413 int v2 = dic3d.Faces[iface * 3 + 2];
414
415 if (v0 >= 0 && static_cast<size_t>(v0) < nCorr &&
416 v1 >= 0 && static_cast<size_t>(v1) < nCorr &&
417 v2 >= 0 && static_cast<size_t>(v2) < nCorr) {
418 face_corr[iface] = std::max({point_corr[v0], point_corr[v1], point_corr[v2]});
419 } else {
420 face_corr[iface] = std::numeric_limits<double>::quiet_NaN();
421 }
422 }
423 }
424 LOG_INFO << " ✓ Face correlation recomputed for " << vertices_all_frames.size() << " frames";
425 } else {
426 LOG_WARN << " ⚠ corrComb not available (empty or frame count mismatch), skipping face correlation";
427 }
428
429 // Compute RBM if enabled
430 std::vector<Utils::RigidTransform> rbm_transforms;
431 std::vector<std::vector<Eigen::Vector3d>> vertices_all_frames_ARBM;
432
433 if (config_.step_f_compute_rbm) {
434 LOG_INFO << "\nComputing rigid body motion (RBM) transformations...";
435 rbm_transforms.resize(vertices_all_frames.size());
436 vertices_all_frames_ARBM.reserve(vertices_all_frames.size());
437
438 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
439 Utils::RigidTransform transform;
440 bool success = Utils::computeRigidTransform(
441 vertices_all_frames[iframe], // from
442 vertices_ref, // to (reference frame)
443 transform
444 );
445
446 rbm_transforms[iframe] = transform;
447
448 if (success) {
449 auto verts_arbm = Utils::applyRigidTransform(vertices_all_frames[iframe], transform);
450 vertices_all_frames_ARBM.push_back(verts_arbm);
451 } else {
452 LOG_WARN << " RBM computation failed for frame " << iframe;
453 vertices_all_frames_ARBM.push_back(vertices_all_frames[iframe]);
454 }
455 }
456 LOG_INFO << " ✓ RBM transformations computed for " << vertices_all_frames.size() << " frames";
457 } else {
458 LOG_INFO << "\nRBM/ARBM computation disabled (step_f_compute_rbm=false)";
459 }
460
461 // Compute 3D deformation and strain
462 LOG_INFO << "\nComputing 3D surface deformation...";
463 LOG_INFO << " Method: Triangular Cosserat Point Elements (TCPE)";
464 LOG_INFO << " Deformation type: Cumulative (reference = frame 1)";
465
467 dic3d.Faces,
468 vertices_ref,
469 vertices_all_frames,
470 true // cumulative: use frame 1 as reference for all frames
471 );
472
473 LOG_INFO << " ✓ Deformation computation complete";
474
475 // Diagnostic: report deformation result quality for first few frames
476 {
477 for (size_t iframe = 0; iframe < std::min(deform_result.n_frames, size_t(3)); ++iframe) {
478 const auto& fr = deform_result.frames[iframe];
479 size_t nan_F = 0, nan_E = 0, nan_J = 0, inf_count = 0, large_count = 0;
480 for (size_t i = 0; i < fr.Fmat.size(); ++i) {
481 if (fr.Fmat[i].array().isNaN().any()) nan_F++;
482 if (fr.Emat.size() > i && fr.Emat[i].array().isNaN().any()) nan_E++;
483 if (fr.J.size() > i && std::isnan(fr.J[i])) nan_J++;
484 if (fr.J.size() > i && std::isinf(fr.J[i])) inf_count++;
485 if (fr.Emgn.size() > i && std::isfinite(fr.Emgn[i]) && std::abs(fr.Emgn[i]) > 1.0)
486 large_count++;
487 }
488 LOG_DEBUG << " [diag] Deform frame " << iframe << ": "
489 << "NaN_F=" << nan_F << " NaN_E=" << nan_E << " NaN_J=" << nan_J
490 << " Inf=" << inf_count << " |Emgn|>1=" << large_count
491 << " / " << fr.Fmat.size() << " faces";
492 }
493 }
494
495 FrameDeformationResult deform_result_ARBM;
496 if (config_.step_f_compute_rbm) {
497 LOG_INFO << "\nComputing 3D surface deformation (after RBM removal)...";
498 deform_result_ARBM = computeTriSurfaceDeformation(
499 dic3d.Faces,
500 vertices_all_frames_ARBM[0],
501 vertices_all_frames_ARBM,
502 true // cumulative
503 );
504 LOG_INFO << " ✓ Deformation computation complete (ARBM)";
505 }
506
507 // Build DIC3DPPresults structure
508 LOG_INFO << "\nBuilding DIC3DPPresults structure...";
509 DIC3DPPresults ppresults;
510
511 // Store full deformation result for MAT file writing (all 36 fields)
512 ppresults.deform_full = deform_result;
513
514 // Compute face isotropy index for each frame
515 LOG_INFO << "\nComputing face isotropy index...";
516 ppresults.FaceIsoInd.resize(vertices_all_frames.size());
517 for (size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
518 ppresults.FaceIsoInd[iframe] = computeFaceIsotropyIndex(
519 dic3d.Faces,
520 vertices_all_frames[iframe]
521 );
522 }
523 LOG_INFO << " ✓ Face isotropy index computed for " << vertices_all_frames.size() << " frames";
524
525 // Copy all fields from DIC3Dcombined (inheritance)
526 ppresults.pairIndices = dic3d.pairIndices;
527 ppresults.Points3D = dic3d.Points3D;
528 ppresults.Faces = dic3d.Faces;
529 ppresults.FaceColors = dic3d.FaceColors;
530 ppresults.corrComb = dic3d.corrComb;
531 ppresults.FaceCorrComb = dic3d.FaceCorrComb;
532 ppresults.FaceCentroids = dic3d.FaceCentroids;
533 ppresults.Disp = dic3d.Disp;
534 ppresults.calibration = dic3d.calibration;
535 ppresults.distortion = dic3d.distortion;
536 ppresults.FacePairInds = dic3d.FacePairInds;
537 ppresults.PointPairInds = dic3d.PointPairInds;
538 ppresults.DIC2Dinfo = dic3d.DIC2Dinfo;
539 ppresults.AllPairsResults = dic3d.AllPairsResults;
540
541 // Set deformation type and frame count
542 ppresults.deftype = "cum"; // cumulative deformation
543 ppresults.n_frames = deform_result.n_frames;
544
545 // Convert FrameDeformationResult to DeformData
546 // This extracts the key fields needed for MATLAB compatibility
547 ppresults.Deform.F.resize(deform_result.n_frames);
548 ppresults.Deform.strain.resize(deform_result.n_frames);
549 ppresults.Deform.princStrain.resize(deform_result.n_frames);
550 ppresults.Deform.maxShearStrain.resize(deform_result.n_frames);
551
552 // nFaces already defined above (line 204)
553
554 for (size_t iframe = 0; iframe < deform_result.n_frames; ++iframe) {
555 const auto& frame = deform_result.frames[iframe];
556
557 // Deformation gradient F (3x3 per face)
558 DeformGradient& F = ppresults.Deform.F[iframe];
559 F.F11.resize(nFaces);
560 F.F12.resize(nFaces);
561 F.F13.resize(nFaces);
562 F.F21.resize(nFaces);
563 F.F22.resize(nFaces);
564 F.F23.resize(nFaces);
565 F.F31.resize(nFaces);
566 F.F32.resize(nFaces);
567 F.F33.resize(nFaces);
568
569 for (size_t iface = 0; iface < nFaces; ++iface) {
570 F.F11[iface] = frame.Fmat[iface](0, 0);
571 F.F12[iface] = frame.Fmat[iface](0, 1);
572 F.F13[iface] = frame.Fmat[iface](0, 2);
573 F.F21[iface] = frame.Fmat[iface](1, 0);
574 F.F22[iface] = frame.Fmat[iface](1, 1);
575 F.F23[iface] = frame.Fmat[iface](1, 2);
576 F.F31[iface] = frame.Fmat[iface](2, 0);
577 F.F32[iface] = frame.Fmat[iface](2, 1);
578 F.F33[iface] = frame.Fmat[iface](2, 2);
579 }
580
581 // Green-Lagrangian strain tensor E (3x3 per face)
582 StrainTensor& E = ppresults.Deform.strain[iframe];
583 E.E11.resize(nFaces);
584 E.E12.resize(nFaces);
585 E.E13.resize(nFaces);
586 E.E21.resize(nFaces);
587 E.E22.resize(nFaces);
588 E.E23.resize(nFaces);
589 E.E31.resize(nFaces);
590 E.E32.resize(nFaces);
591 E.E33.resize(nFaces);
592
593 for (size_t iface = 0; iface < nFaces; ++iface) {
594 E.E11[iface] = frame.Emat[iface](0, 0);
595 E.E12[iface] = frame.Emat[iface](0, 1);
596 E.E13[iface] = frame.Emat[iface](0, 2);
597 E.E21[iface] = frame.Emat[iface](1, 0);
598 E.E22[iface] = frame.Emat[iface](1, 1);
599 E.E23[iface] = frame.Emat[iface](1, 2);
600 E.E31[iface] = frame.Emat[iface](2, 0);
601 E.E32[iface] = frame.Emat[iface](2, 1);
602 E.E33[iface] = frame.Emat[iface](2, 2);
603 }
604
605 // Principal strains (2 per face: Epc1, Epc2)
606 std::vector<double>& princStrain = ppresults.Deform.princStrain[iframe];
607 princStrain.resize(nFaces * 2);
608 for (size_t iface = 0; iface < nFaces; ++iface) {
609 princStrain[iface * 2 + 0] = frame.Epc1[iface];
610 princStrain[iface * 2 + 1] = frame.Epc2[iface];
611 }
612
613 // Max shear strain (1 per face)
614 std::vector<double>& maxShear = ppresults.Deform.maxShearStrain[iframe];
615 maxShear = frame.EShearMax;
616 }
617
618 LOG_INFO << " Populated deformation data for " << ppresults.n_frames << " frames";
619
620 if (config_.step_f_compute_rbm) {
621 // Populate ARBM deformation data (after rigid body motion removal)
622 LOG_INFO << " Populating ARBM deformation data...";
623 ppresults.Deform_ARBM.F.resize(deform_result_ARBM.n_frames);
624 ppresults.Deform_ARBM.strain.resize(deform_result_ARBM.n_frames);
625 ppresults.Deform_ARBM.princStrain.resize(deform_result_ARBM.n_frames);
626 ppresults.Deform_ARBM.maxShearStrain.resize(deform_result_ARBM.n_frames);
627
628 for (size_t iframe = 0; iframe < deform_result_ARBM.n_frames; ++iframe) {
629 const auto& frame = deform_result_ARBM.frames[iframe];
630
631 DeformGradient& F = ppresults.Deform_ARBM.F[iframe];
632 F.F11.resize(nFaces); F.F12.resize(nFaces); F.F13.resize(nFaces);
633 F.F21.resize(nFaces); F.F22.resize(nFaces); F.F23.resize(nFaces);
634 F.F31.resize(nFaces); F.F32.resize(nFaces); F.F33.resize(nFaces);
635
636 for (size_t iface = 0; iface < nFaces; ++iface) {
637 F.F11[iface] = frame.Fmat[iface](0, 0);
638 F.F12[iface] = frame.Fmat[iface](0, 1);
639 F.F13[iface] = frame.Fmat[iface](0, 2);
640 F.F21[iface] = frame.Fmat[iface](1, 0);
641 F.F22[iface] = frame.Fmat[iface](1, 1);
642 F.F23[iface] = frame.Fmat[iface](1, 2);
643 F.F31[iface] = frame.Fmat[iface](2, 0);
644 F.F32[iface] = frame.Fmat[iface](2, 1);
645 F.F33[iface] = frame.Fmat[iface](2, 2);
646 }
647
648 StrainTensor& E = ppresults.Deform_ARBM.strain[iframe];
649 E.E11.resize(nFaces); E.E12.resize(nFaces); E.E13.resize(nFaces);
650 E.E21.resize(nFaces); E.E22.resize(nFaces); E.E23.resize(nFaces);
651 E.E31.resize(nFaces); E.E32.resize(nFaces); E.E33.resize(nFaces);
652
653 for (size_t iface = 0; iface < nFaces; ++iface) {
654 E.E11[iface] = frame.Emat[iface](0, 0);
655 E.E12[iface] = frame.Emat[iface](0, 1);
656 E.E13[iface] = frame.Emat[iface](0, 2);
657 E.E21[iface] = frame.Emat[iface](1, 0);
658 E.E22[iface] = frame.Emat[iface](1, 1);
659 E.E23[iface] = frame.Emat[iface](1, 2);
660 E.E31[iface] = frame.Emat[iface](2, 0);
661 E.E32[iface] = frame.Emat[iface](2, 1);
662 E.E33[iface] = frame.Emat[iface](2, 2);
663 }
664
665 std::vector<double>& princStrain = ppresults.Deform_ARBM.princStrain[iframe];
666 princStrain.resize(nFaces * 2);
667 for (size_t iface = 0; iface < nFaces; ++iface) {
668 princStrain[iface * 2 + 0] = frame.Epc1[iface];
669 princStrain[iface * 2 + 1] = frame.Epc2[iface];
670 }
671 ppresults.Deform_ARBM.maxShearStrain[iframe] = frame.EShearMax;
672 }
673
674 // Populate RBM transformation matrices
675 ppresults.RBM.RotMat.resize(rbm_transforms.size());
676 ppresults.RBM.TransVec.resize(rbm_transforms.size());
677 for (size_t iframe = 0; iframe < rbm_transforms.size(); ++iframe) {
678 const auto& transform = rbm_transforms[iframe];
679 ppresults.RBM.RotMat[iframe].resize(9);
680 for (int i = 0; i < 3; ++i)
681 for (int j = 0; j < 3; ++j)
682 ppresults.RBM.RotMat[iframe][i * 3 + j] = transform.R(i, j);
683 ppresults.RBM.TransVec[iframe] = {transform.t(0), transform.t(1), transform.t(2)};
684 }
685
686 // Populate Points3D_ARBM
687 ppresults.Points3D_ARBM_x.resize(vertices_all_frames_ARBM.size());
688 ppresults.Points3D_ARBM_y.resize(vertices_all_frames_ARBM.size());
689 ppresults.Points3D_ARBM_z.resize(vertices_all_frames_ARBM.size());
690 for (size_t iframe = 0; iframe < vertices_all_frames_ARBM.size(); ++iframe) {
691 const auto& verts = vertices_all_frames_ARBM[iframe];
692 ppresults.Points3D_ARBM_x[iframe].resize(verts.size());
693 ppresults.Points3D_ARBM_y[iframe].resize(verts.size());
694 ppresults.Points3D_ARBM_z[iframe].resize(verts.size());
695 for (size_t i = 0; i < verts.size(); ++i) {
696 ppresults.Points3D_ARBM_x[iframe][i] = verts[i].x();
697 ppresults.Points3D_ARBM_y[iframe][i] = verts[i].y();
698 ppresults.Points3D_ARBM_z[iframe][i] = verts[i].z();
699 }
700 }
701
702 LOG_INFO << " ✓ Populated all deformation data (with RBM and ARBM)";
703 }
704
705 // Save DIC3DPPresults using configured format
706 {
707 auto pp_serializer = cppxdic::DataSerializer::create(config_.data_format);
708 std::string ppout = Utils::buildDic3DPPresultsFilePath(output_dir, config_.num_pair, config_.fileversion, pp_serializer->extension());
709
710 if (std::filesystem::exists(ppout)) {
711 LOG_INFO << "\nCheckpoint found: " << ppout << " (skipping)";
712 } else {
713 LOG_INFO << "\nWriting DIC3DPPresults to: " << ppout;
714 if (!pp_serializer->saveDIC3DPPresults(ppout, ppresults)) {
715 LOG_ERROR << "Failed to write DIC3DPPresults";
716 return false;
717 }
718 LOG_INFO << "✓ Saved DIC3DPPresults: " << ppout;
719 }
720 }
721
722 // Generate visualization exports if enabled
723 if (config_.mapLogic) {
724 LOG_INFO << "\n=== Generating Visualization Exports ===";
725 try {
726 cppxdic::Visualization viz(config_);
727
728 // Print trial info
729 viz.printTrialInfo(ppresults);
730
731 // Apply filtering if enabled
732 viz.applyTemporalFilter(ppresults);
733
734 // Prepare visualization data
735 auto vis_data = viz.prepareVisualizationData(ppresults);
736
737 // Create viz directory if needed
738 std::filesystem::create_directories(output_dir + "viz");
739
740 // Export visualization data
741 auto viz_path = Utils::buildVizPath(output_dir, "trial.txt");
742 viz.exportData(vis_data, viz_path);
743
744 // Generate summary statistics
745 auto stats_path = Utils::buildVizPath(output_dir, "trial_summary.txt");
746 viz.generateSummaryStats(ppresults, stats_path);
747
748 // Generate videos if enabled
749 if (config_.generate_videos) {
750 LOG_INFO << "\n--- Generating Videos ---";
751 for (const auto& field : config_.plotopt) {
752 std::string video_path = output_dir + "viz/" + field + "_video.avi";
753 viz.generateVideo(vis_data, video_path, field);
754 }
755 // Also generate displacement magnitude video if not in plotopt
756 if (std::find(config_.plotopt.begin(), config_.plotopt.end(), "DispMgn") == config_.plotopt.end()) {
757 if (vis_data.FaceScalars.count("DispMgn") || !ppresults.Disp.DispMgn.empty()) {
758 std::string video_path = output_dir + "viz/DispMgn_video.avi";
759 viz.generateVideo(vis_data, video_path, "DispMgn");
760 }
761 }
762 }
763
764 LOG_INFO << "✓ Visualization exports complete";
765
766 } catch (const std::exception& e) {
767 LOG_WARN << "Visualization export failed: " << e.what();
768 LOG_WARN << " (Continuing anyway...)";
769 }
770 }
771
772 LOG_INFO << "\n=== Step F Complete ===";
773 LOG_INFO << "✓ 3D deformation and strain analysis finished for trial " << trial;
774 }
775
776 return true;
777
778 } catch (const std::exception& e) {
779 LOG_ERROR << "in Deformation/Strain Analysis: " << e.what();
780 return false;
781 }
782}
783bool DicAnalysis::dic3DReconstruction(const std::vector<int>& trial_target) {
784 LOG_INFO << "Starting 3D Reconstruction (Step E)...";
785 try {
786
787 auto solve_3d = [](const std::vector<double>& L1, const std::vector<double>& L2, double x1, double y1, double x2, double y2){
788 Eigen::Matrix<double, 4, 3> M;
789 Eigen::Matrix<double, 4, 1> V;
790
791 // Match MATLAB DLT11Reconstruction exactly:
792 // P3D(ii,:) = M \ V
793 M <<
794 x1 * L1[8] - L1[0], x1 * L1[9] - L1[1], x1 * L1[10] - L1[2],
795 y1 * L1[8] - L1[4], y1 * L1[9] - L1[5], y1 * L1[10] - L1[6],
796 x2 * L2[8] - L2[0], x2 * L2[9] - L2[1], x2 * L2[10] - L2[2],
797 y2 * L2[8] - L2[4], y2 * L2[9] - L2[5], y2 * L2[10] - L2[6];
798
799 V <<
800 L1[3] - x1,
801 L1[7] - y1,
802 L2[3] - x2,
803 L2[7] - y2;
804
805 Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 4, 3>> qr(M);
806 if (qr.rank() < 3) {
807 double nan = std::numeric_limits<double>::quiet_NaN();
808 return std::array<double, 3>{nan, nan, nan};
809 }
810
811 Eigen::Matrix<double, 3, 1> X = qr.solve(V);
812 return std::array<double, 3>{X[0], X[1], X[2]};
813 };
814
815 for (int trial : trial_target) {
816 // Collect all individual pair results before stitching
817 std::vector<DIC3DpairResults> all_pairs;
818 std::vector<DIC2DPairResults> dic2d_info;
819
820 // Format trial as 3-digit string (e.g., "005")
821 std::ostringstream trial_str;
822 trial_str << std::setw(3) << std::setfill('0') << trial;
823
824 // ------------------------------------------------------------------
825 // Load DLT calibrations for ALL cameras upfront (matches MATLAB step3)
826 // MATLAB: DLTstructAllCams{ic} = load(DLTpath{ic}).DLTstructCam
827 // ------------------------------------------------------------------
828 std::string calib_dir = Utils::buildCalibDir(config_);
829
830 // Collect unique camera indices across all pairs
831 std::set<int> unique_cams;
832 for (int p = 1; p <= config_.num_pair; ++p) {
833 int c1, c2;
834 Utils::getCamerasForPair(p, c1, c2);
835 unique_cams.insert(c1);
836 unique_cams.insert(c2);
837 }
838
839 // Load DLT calibration for each unique camera
840 std::map<int, DLTCalibrationData> dlt_all_cams;
841 for (int cam_id : unique_cams) {
842 // Match MATLAB pattern: DLTstruct_cam_<id>.mat (or *cam_<id>.mat via glob)
843 std::string dlt_path;
844 if (std::filesystem::exists(calib_dir)) {
845 std::string search_suffix = "cam_" + std::to_string(cam_id) + ".mat";
846 for (const auto& entry : std::filesystem::directory_iterator(calib_dir)) {
847 if (!entry.is_regular_file()) continue;
848 auto name = entry.path().filename().string();
849 if (name.size() >= search_suffix.size() &&
850 name.compare(name.size() - search_suffix.size(), search_suffix.size(), search_suffix) == 0) {
851 dlt_path = entry.path().string();
852 break;
853 }
854 }
855 }
856
857 if (dlt_path.empty()) {
858 LOG_ERROR << "DLT calibration file not found for camera " << cam_id
859 << " in " << calib_dir;
860 continue;
861 }
862
863 DLTCalibrationData calib;
864 if (cppxdic::MatReader::loadDLTCalibration(dlt_path, calib)) {
865 dlt_all_cams[cam_id] = std::move(calib);
866 } else {
867 LOG_ERROR << "Failed to load DLT calibration for camera " << cam_id;
868 }
869 }
870
871 if (dlt_all_cams.empty()) {
872 LOG_ERROR << "No DLT calibrations loaded. Cannot proceed with 3D reconstruction.";
873 continue;
874 }
875 LOG_INFO << "Loaded DLT calibrations for " << dlt_all_cams.size() << " cameras";
876
877 // ------------------------------------------------------------------
878 // Load distortion parameters for all cameras (if available)
879 // ------------------------------------------------------------------
880 std::map<int, Utils::CameraParameters> distortion_all_cams;
881 std::map<int, std::string> distortion_paths;
882 if (std::filesystem::exists(calib_dir)) {
883 for (int cam_id : unique_cams) {
884 std::string search_str = "cam_" + std::to_string(cam_id);
885 for (const auto& entry : std::filesystem::directory_iterator(calib_dir)) {
886 if (!entry.is_regular_file()) continue;
887 auto name = entry.path().filename().string();
888 std::string lower = name;
889 std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower);
890 if (lower.find("cameracbparameters") != std::string::npos &&
891 lower.find(search_str) != std::string::npos) {
893 if (Utils::loadCameraParameters(entry.path().string(), params)) {
894 distortion_all_cams[cam_id] = params;
895 distortion_paths[cam_id] = entry.path().string();
896 }
897 break;
898 }
899 }
900 }
901 }
902 bool use_distortion_removal = config_.step_e_distortion_removal;
903 if (use_distortion_removal) {
904 LOG_INFO << "Distortion removal enabled for all " << unique_cams.size() << " cameras";
905
906 // ------------------------------------------------------------------
907 // Recalculate DLT parameters from undistorted calibration centroids
908 // Matches MATLAB step3 lines 84-98:
909 // P2Dtemp = DLTstructAllCams{ic}.imageCentroids;
910 // [P2Dtemp] = undistortPoints(P2Dtemp, distortionPar{ic});
911 // C3Dtemp = reshape(C3Dtrue(:,columns,:), nRows*nCols, 3);
912 // L = DLT11Calibration(P2Dtemp, C3Dtemp);
913 // DLTstructAllCams{ic}.DLTparams = L;
914 // ------------------------------------------------------------------
915 for (int cam_id : unique_cams) {
916 auto& dlt = dlt_all_cams[cam_id];
917 const auto& dist = distortion_all_cams[cam_id];
918
919 if (dlt.imageCentroids.empty() || dlt.C3Dtrue.empty() || dlt.columns.empty()) {
920 LOG_ERROR << "Camera " << cam_id << ": missing calibration data for DLT recalculation, "
921 << "using original DLT params with per-point undistortion fallback";
922 continue;
923 }
924
925 // Undistort calibration image centroids (Nx2 row-major)
926 size_t nCentroids = dlt.imageCentroids_rows;
927 std::vector<cv::Point2d> centroids_in(nCentroids);
928 for (size_t i = 0; i < nCentroids; ++i) {
929 centroids_in[i] = cv::Point2d(dlt.imageCentroids[i * 2 + 0],
930 dlt.imageCentroids[i * 2 + 1]);
931 }
932 std::vector<cv::Point2d> centroids_undist;
933 Utils::undistortPoints(centroids_in, dist, centroids_undist);
934
935 // Convert undistorted centroids back to flat array (Nx2)
936 std::vector<double> P2D(nCentroids * 2);
937 for (size_t i = 0; i < nCentroids; ++i) {
938 P2D[i * 2 + 0] = centroids_undist[i].x;
939 P2D[i * 2 + 1] = centroids_undist[i].y;
940 }
941
942 // Reshape C3Dtrue(:, columns, :) -> (nRows * nSelectedCols) x 3
943 // C3Dtrue is stored row-major as (dim0 x dim1 x 3)
944 // columns contains 1-based column indices used in calibration
945 size_t d0 = dlt.C3Dtrue_dim0; // nRows
946 size_t nCols = dlt.columns.size();
947 std::vector<double> P3D(d0 * nCols * 3);
948 for (size_t ci = 0; ci < nCols; ++ci) {
949 size_t col = static_cast<size_t>(dlt.columns[ci]) - 1; // MATLAB 1-based to 0-based
950 for (size_t row = 0; row < d0; ++row) {
951 size_t dst_idx = (ci * d0 + row); // column-major order matching MATLAB reshape
952 size_t src_idx = (row * dlt.C3Dtrue_dim1 + col); // row-major C3Dtrue indexing
953 P3D[dst_idx * 3 + 0] = dlt.C3Dtrue[src_idx * 3 + 0];
954 P3D[dst_idx * 3 + 1] = dlt.C3Dtrue[src_idx * 3 + 1];
955 P3D[dst_idx * 3 + 2] = dlt.C3Dtrue[src_idx * 3 + 2];
956 }
957 }
958
959 size_t nPoints = d0 * nCols;
960 if (nPoints != nCentroids) {
961 LOG_ERROR << "Camera " << cam_id << ": centroid count (" << nCentroids
962 << ") != C3D point count (" << nPoints << "), skipping DLT recalculation";
963 continue;
964 }
965
966 // Recalculate DLT parameters
967 std::vector<double> L_new;
968 if (Utils::DLT11Calibration(P2D.data(), P3D.data(), nPoints, L_new)) {
969 dlt.DLTparams = L_new;
970 LOG_INFO << "Camera " << cam_id << ": DLT params recalculated from "
971 << nPoints << " undistorted calibration points";
972 } else {
973 LOG_ERROR << "Camera " << cam_id << ": DLT11Calibration failed";
974 }
975 }
976 } else if (!distortion_all_cams.empty()) {
977 LOG_INFO << "Distortion parameters found for " << distortion_all_cams.size()
978 << "/" << unique_cams.size() << " cameras (disabled - need all)";
979 use_distortion_removal = false;
980 }
981
982 // ------------------------------------------------------------------
983 // Process each stereo pair
984 // ------------------------------------------------------------------
985 for (int pair = 1; pair <= config_.num_pair; ++pair) {
986 LOG_INFO << "\n=== Processing Pair " << pair << " ===";
987
988 // Get camera numbers for this pair
989 int cam_1, cam_2;
990 Utils::getCamerasForPair(pair, cam_1, cam_2);
991
992 // Look up DLT parameters for this pair's cameras
993 if (dlt_all_cams.find(cam_1) == dlt_all_cams.end() ||
994 dlt_all_cams.find(cam_2) == dlt_all_cams.end()) {
995 LOG_ERROR << "Missing DLT calibration for pair " << pair
996 << " (cam " << cam_1 << " or " << cam_2 << "). Skipping.";
997 continue;
998 }
999 const auto& dlt_cam1 = dlt_all_cams[cam_1];
1000 const auto& dlt_cam2 = dlt_all_cams[cam_2];
1001 const std::vector<double>& L1 = dlt_cam1.DLTparams;
1002 const std::vector<double>& L2 = dlt_cam2.DLTparams;
1003
1004 // ---------------------------------------------------------------
1005 // Load DIC2DpairResults (MATLAB step3 lines 39-47)
1006 // All 2D point extraction, matching, faces, colors, and correlation
1007 // are already computed by formatOutput (step2_dic_finish equivalent).
1008 // ---------------------------------------------------------------
1009 auto d_serializer = cppxdic::DataSerializer::create(config_.data_format);
1010 std::string output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
1011 std::string dic2d_file = Utils::buildDic2DPairResultsFilePath(output_dir, cam_1, cam_2, d_serializer->extension());
1012
1013 if (!std::filesystem::exists(dic2d_file)) {
1014 LOG_ERROR << "Missing DIC2DPairResults for trial " << trial << ", pair " << pair << ". Skipping.";
1015 LOG_ERROR << " Expected: " << dic2d_file;
1016 continue;
1017 }
1018
1019 DIC2DPairResults DIC2D;
1020 if (!d_serializer->loadDIC2DPairResults(dic2d_file, DIC2D)) {
1021 LOG_ERROR << "Failed to load DIC2DPairResults: " << dic2d_file << ". Skipping.";
1022 continue;
1023 }
1024
1025 if (config_.debug_mode) {
1026 LOG_DEBUG << "Loaded DIC2DPairResults from: " << dic2d_file;
1027 LOG_DEBUG << " nImages=" << DIC2D.nImages << ", Points=" << DIC2D.Points.size()
1028 << ", Faces=" << DIC2D.Faces.size()/3;
1029 LOG_DEBUG << " DLT cam" << cam_1 << ": " << dlt_cam1.filePath;
1030 LOG_DEBUG << " DLT cam" << cam_2 << ": " << dlt_cam2.filePath;
1031 }
1032
1033 // Extract information from 2D-DIC results (MATLAB step3 lines 130-134)
1034 int nImages = DIC2D.nImages;
1035 const auto& CorCoeff = DIC2D.CorCoeffVec;
1036 const auto& F = DIC2D.Faces;
1037 const auto& FC = DIC2D.FaceColors;
1038 const auto& Points = DIC2D.Points;
1039
1040 // Create individual pair result structure (MATLAB step3 lines 110-136)
1041 DIC3DpairResults pair_result;
1042 pair_result.cameraPairInd = {cam_1, cam_2};
1043 pair_result.DLTpath = {dlt_cam1.filePath, dlt_cam2.filePath};
1044 pair_result.DLTparameters = {L1, L2};
1045 pair_result.Faces = F;
1046 pair_result.FaceColors = FC;
1047
1048 // Set distortion info (MATLAB step3 lines 138-186)
1049 if (use_distortion_removal) {
1050 pair_result.distortionModel = {"distortion", "distortion"};
1051 pair_result.distortionPath = {distortion_paths[cam_1], distortion_paths[cam_2]};
1052 } else {
1053 pair_result.distortionModel = {"none", "none"};
1054 pair_result.distortionPath = {"none", "none"};
1055 }
1056
1057 size_t nFaces = F.size() / 3;
1058 std::vector<double> P3D_ref; // frame 1 reference for displacement
1059
1060 // Loop over images/frames (MATLAB step3 lines 196-221)
1061 for (int ii = 0; ii < nImages; ++ii) {
1062 // Correlated points from 2 cameras (MATLAB: P1=Points{ii}, P2=Points{ii+nImages})
1063 const auto& P1 = Points[ii];
1064 const auto& P2 = Points[ii + nImages];
1065 size_t N = P1.x.size();
1066
1067 // Solve the DLT system (MATLAB: P3D = DLT11Reconstruction(P1, P2, L1, L2))
1068 std::vector<double> pts3d(N * 3, std::numeric_limits<double>::quiet_NaN());
1069 Point3DFrame frame_pts;
1070 frame_pts.x.resize(N); frame_pts.y.resize(N); frame_pts.z.resize(N);
1071 for (size_t k = 0; k < N; ++k) {
1072 double x1 = P1.x[k], y1 = P1.y[k];
1073 double x2 = P2.x[k], y2 = P2.y[k];
1074 if (!std::isnan(x1) && !std::isnan(y1) && !std::isnan(x2) && !std::isnan(y2)) {
1075 auto X = solve_3d(L1, L2, x1, y1, x2, y2);
1076 pts3d[k*3+0] = X[0]; pts3d[k*3+1] = X[1]; pts3d[k*3+2] = X[2];
1077 frame_pts.x[k] = X[0]; frame_pts.y[k] = X[1]; frame_pts.z[k] = X[2];
1078 } else {
1079 frame_pts.x[k] = std::numeric_limits<double>::quiet_NaN();
1080 frame_pts.y[k] = std::numeric_limits<double>::quiet_NaN();
1081 frame_pts.z[k] = std::numeric_limits<double>::quiet_NaN();
1082 }
1083 }
1084 pair_result.Points3D.push_back(frame_pts);
1085
1086 // Combined correlation coefficients (MATLAB: corrComb = max([CorCoeff{ii} CorCoeff{ii+nImages}], [], 2))
1087 std::vector<double> corr_comb(N, 0.0);
1088 if (ii < (int)CorCoeff.size() && (ii + nImages) < (int)CorCoeff.size()) {
1089 for (size_t k = 0; k < N; ++k) {
1090 double cc1 = (k < CorCoeff[ii].size()) ? CorCoeff[ii][k] : 0.0;
1091 double cc2 = (k < CorCoeff[ii + nImages].size()) ? CorCoeff[ii + nImages][k] : 0.0;
1092 corr_comb[k] = std::max(cc1, cc2);
1093 }
1094 }
1095 pair_result.corrComb.push_back(corr_comb);
1096
1097 // Face correlation coefficient (MATLAB: FaceCorrComb = max(corrComb(F), [], 2))
1098 std::vector<double> face_corr(nFaces, std::numeric_limits<double>::quiet_NaN());
1099 for (size_t iface = 0; iface < nFaces; ++iface) {
1100 int v0 = F[iface*3], v1 = F[iface*3+1], v2 = F[iface*3+2];
1101 if (v0 >= 0 && v0 < (int)N && v1 >= 0 && v1 < (int)N && v2 >= 0 && v2 < (int)N) {
1102 face_corr[iface] = std::max({corr_comb[v0], corr_comb[v1], corr_comb[v2]});
1103 }
1104 }
1105 pair_result.FaceCorrComb.push_back(face_corr);
1106
1107 // Compute face centroids (MATLAB: FaceCentroids(iface,:) = mean(P3D(F(iface,:),:)))
1108 std::vector<double> face_centroids(nFaces * 3, std::numeric_limits<double>::quiet_NaN());
1109 for (size_t iface = 0; iface < nFaces; ++iface) {
1110 int v0 = F[iface*3], v1 = F[iface*3+1], v2 = F[iface*3+2];
1111 if (v0 >= 0 && v0 < (int)N && v1 >= 0 && v1 < (int)N && v2 >= 0 && v2 < (int)N) {
1112 face_centroids[iface*3+0] = (pts3d[v0*3+0] + pts3d[v1*3+0] + pts3d[v2*3+0]) / 3.0;
1113 face_centroids[iface*3+1] = (pts3d[v0*3+1] + pts3d[v1*3+1] + pts3d[v2*3+1]) / 3.0;
1114 face_centroids[iface*3+2] = (pts3d[v0*3+2] + pts3d[v1*3+2] + pts3d[v2*3+2]) / 3.0;
1115 }
1116 }
1117 pair_result.FaceCentroids.push_back(face_centroids);
1118
1119 // Compute displacements between frames (MATLAB: DispVec = Points3D{ii} - Points3D{1})
1120 if (ii == 0) P3D_ref = pts3d;
1121 std::vector<double> dispvec(N * 3, std::numeric_limits<double>::quiet_NaN());
1122 std::vector<double> dispmgn(N, std::numeric_limits<double>::quiet_NaN());
1123 for (size_t k = 0; k < N; ++k) {
1124 double dx = pts3d[k*3+0] - P3D_ref[k*3+0];
1125 double dy = pts3d[k*3+1] - P3D_ref[k*3+1];
1126 double dz = pts3d[k*3+2] - P3D_ref[k*3+2];
1127 dispvec[k*3+0] = dx; dispvec[k*3+1] = dy; dispvec[k*3+2] = dz;
1128 dispmgn[k] = std::sqrt(dx*dx + dy*dy + dz*dz);
1129 }
1130 pair_result.Disp.DispVec.push_back(dispvec);
1131 pair_result.Disp.DispMgn.push_back(dispmgn);
1132 }
1133
1134 // Store this pair's result and DIC2D data
1135 all_pairs.push_back(pair_result);
1136 dic2d_info.push_back(std::move(DIC2D));
1137
1138 LOG_INFO << "✓ Pair " << pair << " complete: "
1139 << pair_result.Points3D[0].x.size() << " points, "
1140 << nFaces << " faces";
1141 }
1142
1143 // Stitch all pairs together
1144 LOG_INFO << "\n=== Stitching " << all_pairs.size() << " pairs ===";
1145 DIC3Dcombined stitched;
1146 if (all_pairs.empty()) {
1147 LOG_ERROR << "No pairs successfully reconstructed for trial " << trial;
1148 continue;
1149 } else {
1150 std::vector<int> stitch_pair_order;
1151 bool stitch_pair_forced = false;
1152 bool have_stitch_metadata = false;
1153 bool stitch_metadata_mismatch = false;
1154
1155 for (const auto& dic2d : dic2d_info) {
1156 if (dic2d.pairOrder.empty()) {
1157 LOG_ERROR << "Empty pairOrder in DIC2D for trial " << trial;
1158 continue;
1159 }
1160 if (!have_stitch_metadata) {
1161 LOG_WARN << "Using first pair's stitch metadata for trial " << trial;
1162 stitch_pair_order = dic2d.pairOrder;
1163 stitch_pair_forced = dic2d.pairForced;
1164 have_stitch_metadata = true;
1165 } else if (dic2d.pairOrder != stitch_pair_order ||
1166 dic2d.pairForced != stitch_pair_forced) {
1167 LOG_ERROR << "Inconsistent DIC2D stitch metadata across pairs for trial "
1168 << trial << ". Skipping trial.";
1169 stitch_metadata_mismatch = true;
1170 break;
1171 }
1172 }
1173
1174 if (stitch_metadata_mismatch) {
1175 LOG_WARN << "Skipping trial " << trial << " due to stitch metadata mismatch";
1176 continue;
1177 }
1178
1180 if (all_pairs.size() > 1 && config_.step_e_stitch_mode == "simple") {
1181 LOG_INFO << " step_e_stitch_mode=simple: appending " << all_pairs.size()
1182 << " pairs without overlap removal";
1183 stitched = stitchPairsSimple(all_pairs);
1184 } else if (all_pairs.size() > 1 && all_pairs.size() == static_cast<size_t>(config_.num_pair)) {
1185 if (!have_stitch_metadata) {
1186 bool protocol_loaded = false;
1187 bool protocol_available = false;
1188 cppxdic::ProtocolFileData protocol_data;
1189 stitch_pair_order.reserve(config_.num_pair);
1190 for (int pair = 1; pair <= config_.num_pair; ++pair) {
1191 stitch_pair_order.push_back(pair);
1192 }
1193
1194 try {
1195 std::string protocol_dir = Utils::buildProtocolDir(config_, true, true, true, true);
1196 auto protocol_files = Utils::findFiles(protocol_dir, "*.mat");
1197 if (!protocol_files.empty()) {
1198 protocol_loaded = true;
1199 protocol_available = MatReader::loadProtocol(protocol_files.front(), protocol_data);
1200 }
1201 } catch (...) {
1202 protocol_available = false;
1203 }
1204
1205 if (protocol_loaded && protocol_available &&
1206 trial > 0 && static_cast<size_t>(trial) <= protocol_data.trials.size()) {
1207 const std::string& direction = protocol_data.trials[trial - 1].direction;
1208 if (direction == "Ubnf") {
1209 stitch_pair_order = {2, 1};
1210 stitch_pair_forced = true;
1211 } else if (direction == "Rbnf") {
1212 stitch_pair_order = {1, 2};
1213 stitch_pair_forced = true;
1214 } else {
1215 stitch_pair_forced = false;
1216 }
1217 for (int pair = 1; pair <= config_.num_pair; ++pair) {
1218 if (std::find(stitch_pair_order.begin(), stitch_pair_order.end(), pair) == stitch_pair_order.end()) {
1219 stitch_pair_order.push_back(pair);
1220 }
1221 }
1222 }
1223 LOG_INFO << " Stitch metadata missing from DIC2D results; using legacy fallback";
1224 }
1225
1226 {
1227 std::ostringstream stitch_order_oss;
1228 for (size_t i = 0; i < stitch_pair_order.size(); ++i) {
1229 if (i > 0) stitch_order_oss << ", ";
1230 stitch_order_oss << stitch_pair_order[i];
1231 }
1232 LOG_INFO << " Stitch order: " << stitch_order_oss.str();
1233 }
1234 LOG_INFO << " Pair forced metadata: " << (stitch_pair_forced ? "true" : "false");
1235 LOG_INFO << " Ray-cast prefilter (step_e_ray_bvh): "
1236 << (config_.step_e_ray_bvh ? "on (AABB hierarchy)" : "off (full scan)");
1237 stitched = stitchPairsGeometric(all_pairs, stitch_pair_order, stitch_pair_forced);
1238 LOG_INFO << "Geometric Stitching done!";
1239 } else {
1240 if (all_pairs.size() > 1) {
1241 LOG_INFO << " Falling back to simple append stitching because only "
1242 << all_pairs.size() << "/" << config_.num_pair
1243 << " pairs were reconstructed";
1244 }
1245 stitched = stitchPairsSimple(all_pairs);
1246 }
1247
1248 // Store individual pair results in stitched structure
1249 stitched.AllPairsResults = all_pairs;
1250
1251 // DIC2D pair results already collected during per-pair processing (MATLAB step3 line 268)
1252 stitched.DIC2Dinfo = dic2d_info;
1253 LOG_INFO << " Stored " << stitched.DIC2Dinfo.size() << " DIC2D pair results";
1254 }
1255
1256 // Save DIC3Dcombined using configured format
1257 {
1258 auto e_serializer = cppxdic::DataSerializer::create(config_.data_format);
1259 auto output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
1260 std::string e_out = Utils::buildDic3DCombinedFilePath(output_dir, config_.num_pair, e_serializer->extension());
1261
1262 if (std::filesystem::exists(e_out)) {
1263 LOG_INFO << "Checkpoint found: " << e_out << " (skipping)";
1264 } else {
1265 LOG_INFO << "Writing DIC3Dcombined to: " << e_out;
1266 // Crash-safe checkpoint: write to a .partial file and
1267 // rename on success, so a crash mid-save can never leave a
1268 // truncated file that later runs mistake for a valid
1269 // checkpoint (the existence check above cannot tell).
1270 const std::string e_tmp = e_out + ".partial";
1271 if (!e_serializer->saveDIC3Dcombined(e_tmp, stitched)) {
1272 LOG_ERROR << "Failed to write DIC3Dcombined";
1273 std::error_code ec;
1274 std::filesystem::remove(e_tmp, ec);
1275 } else {
1276 std::error_code ec;
1277 std::filesystem::rename(e_tmp, e_out, ec);
1278 if (ec) {
1279 LOG_ERROR << "Failed to finalize DIC3Dcombined: " << ec.message();
1280 } else {
1281 LOG_INFO << "✓ Saved DIC3Dcombined: " << e_out;
1282 }
1283 }
1284 }
1285 }
1286 }
1287 return true;
1288 } catch (const std::exception& e) {
1289 LOG_ERROR << "in 3D Reconstruction: " << e.what();
1290 return false;
1291 }
1292}
1293
1294bool DicAnalysis::setupNcorrAnalysis(const std::vector<std::string>& images,
1295 const std::string& roi_mask_image_path,
1296 DIC_analysis_input& dic_input) {
1297 try {
1298 if (images.size() < 2) {
1299 LOG_ERROR << "Need at least 2 images for DIC analysis";
1300 return false;
1301 }
1302
1303 std::vector<Image2D> ncorr_images;
1304 for (const auto& img_path : images) {
1305 ncorr_images.emplace_back(img_path);
1306 }
1307
1308 // Build ROI from provided mask image (non-zero pixels inside ROI)
1309 Image2D roi_mask(roi_mask_image_path);
1310 ROI2D roi(roi_mask.get_gs() > 0.5);
1311
1312 dic_input = DIC_analysis_input(
1313 ncorr_images,
1314 roi,
1315 3,
1316 INTERP::QUINTIC_BSPLINE_PRECOMPUTE,
1317 SUBREGION::CIRCLE,
1318 config_.subregion_radius,
1319 4,
1320 DIC_analysis_config::NO_UPDATE,
1321 config_.debug_mode
1322 );
1323
1324 return true;
1325
1326 } catch (const std::exception& e) {
1327 LOG_ERROR << "setting up ncorr analysis with ROI mask: " << e.what();
1328 return false;
1329 }
1330}
1331
1333 // Default trial target (equivalent to search_trial2target).
1334 // Trial-level drivers should call run(trial_target) with an explicit list.
1335 std::vector<int> trial_target = {7}; //{7, 12, 25};//#searchTrialTarget();
1336 return run(trial_target);
1337}
1338
1339bool DicAnalysis::run(const std::vector<int>& trial_target) {
1340 // Default: full pipeline (all stages). Preserves historical behaviour.
1341 return run(trial_target, StagePlan{});
1342}
1343
1344bool DicAnalysis::run(const std::vector<int>& trial_target, const StagePlan& plan) {
1345 {
1346 std::ostringstream trial_target_oss;
1347 trial_target_oss << "[";
1348 for (size_t i = 0; i < trial_target.size(); ++i) {
1349 trial_target_oss << trial_target[i];
1350 if (i < trial_target.size() - 1) trial_target_oss << ", ";
1351 }
1352 trial_target_oss << "]";
1353 LOG_INFO << "Trial target set: " << trial_target_oss.str();
1354 }
1355
1356 // Create serializer for format-aware checkpoint checks
1357 auto serializer = cppxdic::DataSerializer::create(config_.data_format);
1358 std::string ext = serializer->extension();
1359
1360 // Helper lambda: Check if 2D DIC outputs exist for all trials/pairs
1361 auto check_2d_outputs_exist = [&]() -> bool {
1362 for (int trial : trial_target) {
1363 auto output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
1364 for (int pair = 1; pair <= config_.num_pair; ++pair) {
1365 int cam1, cam2;
1366 Utils::getCamerasForPair(pair, cam1, cam2);
1367
1368 std::string path = Utils::buildDic2DPairResultsFilePath(output_dir, cam1, cam2, ext);
1369 if (!std::filesystem::exists(path)) {
1370 return false;
1371 }
1372 }
1373 }
1374 return true;
1375 };
1376
1377 // Helper lambda: Check if 3D reconstruction outputs exist
1378 auto check_3d_outputs_exist = [&]() -> bool {
1379 for (int trial : trial_target) {
1380 auto output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
1381 std::string dic3d_file = Utils::buildDic3DCombinedFilePath(output_dir, config_.num_pair, ext);
1382 if (!std::filesystem::exists(dic3d_file)) {
1383 return false;
1384 }
1385 }
1386 return true;
1387 };
1388
1389 // Helper lambda: Check if deformation analysis outputs exist
1390 auto check_deformation_outputs_exist = [&]() -> bool {
1391 for (int trial : trial_target) {
1392 auto output_dir = Utils::buildOutputUntilPhaseDir(config_, trial);
1393 std::string pp_file = Utils::buildDic3DPPresultsFilePath(output_dir, config_.num_pair, config_.fileversion, ext);
1394 if (!std::filesystem::exists(pp_file)) {
1395 return false;
1396 }
1397 }
1398 return true;
1399 };
1400
1401 // STEP D: 2D-DIC (matching / tracking / format) — only if the plan asks for it.
1402 if (plan.anyStepD()) {
1403 bool step_d_complete = check_2d_outputs_exist();
1404 if (step_d_complete) {
1405 LOG_INFO << "\n=== STEP D: 2D-DIC ===";
1406 LOG_INFO << "✓ Checkpoint detected: All 2D DIC output files exist";
1407 LOG_INFO << " Skipping 2D analysis (use existing results)";
1408 } else {
1409 LOG_INFO << "\n=== STEP D: 2D-DIC ===";
1410 LOG_INFO << "Running 2D DIC analysis (stages:"
1411 << (plan.match ? " match" : "") << (plan.track ? " track" : "")
1412 << (plan.format ? " format" : "") << ")...";
1413 auto start = std::chrono::high_resolution_clock::now();
1414 bool success = dic2DAnalysis(trial_target, plan);
1415 auto end = std::chrono::high_resolution_clock::now();
1416 auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
1417
1418 if (!success) {
1419 LOG_ERROR << "2D DIC Analysis failed!";
1420 return false;
1421 }
1422
1423 LOG_INFO << "✓ DIC 2D Analysis done in " << duration.count() / 1000.0 << " s";
1424 }
1425 }
1426
1427 // STEP E: 3D Reconstruction
1428 if (plan.recon) {
1429 bool step_e_complete = check_3d_outputs_exist();
1430 if (step_e_complete) {
1431 LOG_INFO << "\n=== STEP E: 3D Reconstruction ===";
1432 LOG_INFO << "✓ Checkpoint detected: All 3D reconstruction output files exist";
1433 LOG_INFO << " Skipping 3D reconstruction (use existing results)";
1434 } else {
1435 LOG_INFO << "\n=== STEP E: 3D Reconstruction ===";
1436 LOG_INFO << "Running 3D reconstruction...";
1437 auto start = std::chrono::high_resolution_clock::now();
1438 bool success = dic3DReconstruction(trial_target);
1439 auto end = std::chrono::high_resolution_clock::now();
1440 auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
1441
1442 if (!success) {
1443 LOG_ERROR << "3D Reconstruction failed!";
1444 return false;
1445 }
1446
1447 LOG_INFO << "✓ DIC 3D Reconstruction done in " << duration.count() / 1000.0 << " s";
1448 }
1449 }
1450
1451 // STEP F: Deformation analysis
1452 if (plan.deform) {
1453 bool step_f_complete = check_deformation_outputs_exist();
1454 if (step_f_complete) {
1455 LOG_INFO << "\n=== STEP F: Deformation Analysis ===";
1456 LOG_INFO << "✓ Checkpoint detected: Deformation analysis output exists";
1457 LOG_INFO << " Skipping deformation analysis (use existing results)";
1458 } else {
1459 LOG_INFO << "\n=== STEP F: Deformation Analysis ===";
1460 LOG_INFO << "Running deformation analysis...";
1461 auto start = std::chrono::high_resolution_clock::now();
1462 bool success = dicDeformationAnalysis(trial_target);
1463 auto end = std::chrono::high_resolution_clock::now();
1464 auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
1465
1466 if (!success) {
1467 LOG_ERROR << "Deformation Analysis failed!";
1468 return false;
1469 }
1470
1471 LOG_INFO << "✓ DIC Deformation Analysis done in " << duration.count() / 1000.0 << " s";
1472 }
1473 }
1474
1475 LOG_INFO << "\n========================================";
1476 LOG_INFO << "✓ All DIC analysis steps complete!";
1477 LOG_INFO << "========================================";
1478
1479 return true;
1480}
1481
1482std::vector<int> DicAnalysis::searchTrialTarget() {
1483 // Preserve original behaviour: look up trials for the configured subject.
1484 return searchTrialTarget(config_.subject_id);
1485}
1486
1487std::vector<int> DicAnalysis::searchTrialTarget(const std::string& subject) {
1488 std::vector<int> trials;
1489
1490 // Work on a copy of the configuration with the requested subject so that
1491 // arbitrary subjects can be queried without mutating config_.
1492 Config cfg = config_;
1493 cfg.overrideSubject(subject);
1494 cfg.updateVariables();
1495
1496 // Build protocol directory path
1497 std::string protocol_dir = Utils::buildProtocolDir(cfg, true, true, true, true);
1498
1499 // Find protocol .mat file
1500 auto protos = Utils::findFiles(protocol_dir, "*.mat");
1501 if (protos.empty()) {
1502 LOG_ERROR << "Protocol file not found in: " << protocol_dir;
1503 // Fallback to reference + next
1504 trials.push_back(cfg.ref_trial_id);
1505 trials.push_back(cfg.ref_trial_id + 1);
1506 return trials;
1507 }
1508
1509 std::string proto_file = protos.front();
1510
1511 // Open MAT file
1512 mat_t *matfp = Mat_Open(proto_file.c_str(), MAT_ACC_RDONLY);
1513 if (!matfp) {
1514 LOG_ERROR << "Failed to open MAT file: " << proto_file;
1515 trials.push_back(cfg.ref_trial_id);
1516 trials.push_back(cfg.ref_trial_id + 1);
1517 return trials;
1518 }
1519
1520 // Load 'cond' struct
1521 matvar_t *cond = Mat_VarRead(matfp, "cond");
1522 if (!cond || cond->class_type != MAT_C_STRUCT) {
1523 if (cond) Mat_VarFree(cond);
1524 Mat_Close(matfp);
1525 LOG_ERROR << "Variable 'cond' not found or not a struct in: " << proto_file;
1526 trials.push_back(cfg.ref_trial_id);
1527 trials.push_back(cfg.ref_trial_id + 1);
1528 return trials;
1529 }
1530
1531 // Read titles (cell array of strings)
1532 matvar_t *titles = Mat_VarGetStructFieldByName(cond, "titles", 0);
1533 matvar_t *table = Mat_VarGetStructFieldByName(cond, "table", 0);
1534 if (!titles || titles->class_type != MAT_C_CELL || !table || table->class_type != MAT_C_CELL) {
1535 if (cond) Mat_VarFree(cond);
1536 Mat_Close(matfp);
1537 LOG_ERROR << "Fields 'titles' or 'table' missing or of wrong type in 'cond'";
1538 trials.push_back(cfg.ref_trial_id);
1539 trials.push_back(cfg.ref_trial_id + 1);
1540 return trials;
1541 }
1542
1543 // Extract column names
1544 std::vector<std::string> col_names;
1545 size_t ncols = titles->dims[1];
1546 col_names.reserve(ncols);
1547 for (size_t j = 0; j < ncols; ++j) {
1548 matvar_t *cell = static_cast<matvar_t **>(titles->data)[j];
1549 std::string name;
1550 if (cell && cell->class_type == MAT_C_CHAR && cell->data) {
1551 size_t len = cell->nbytes / cell->data_size;
1552 name.assign(static_cast<const char *>(cell->data), len);
1553 // Titles might have trailing nulls; trim
1554 while (!name.empty() && (name.back() == '\0' || name.back() == ' ')) name.pop_back();
1555 }
1556 col_names.push_back(name);
1557 }
1558
1559 // Identify indices for 'dir', 'nf', 'spddxl'
1560 auto find_col = [&](const std::string &key) -> int {
1561 for (size_t j = 0; j < col_names.size(); ++j) {
1562 if (col_names[j] == key) return static_cast<int>(j);
1563 }
1564 return -1;
1565 };
1566 int dir_idx = find_col("dir");
1567 int nf_idx = find_col("nf");
1568 int spd_idx = find_col("spddxl");
1569 if (dir_idx < 0 || nf_idx < 0 || spd_idx < 0) {
1570 if (cond) Mat_VarFree(cond);
1571 Mat_Close(matfp);
1572 LOG_ERROR << "Required columns not found in titles (need 'dir','nf','spddxl')";
1573 trials.push_back(cfg.ref_trial_id);
1574 trials.push_back(cfg.ref_trial_id + 1);
1575 return trials;
1576 }
1577
1578 // Table dimensions: Ntrial x Ncond
1579 size_t ntrial = table->dims[0];
1580 size_t ncond = table->dims[1];
1581
1582 // Pre-extract columns from cell table
1583 auto cell_at = [&](size_t i, size_t j) -> matvar_t * {
1584 size_t idx = i + j * ntrial; // column-major
1585 return static_cast<matvar_t **>(table->data)[idx];
1586 };
1587
1588 std::vector<std::string> dircol(ntrial);
1589 std::vector<double> nfcol(ntrial, std::numeric_limits<double>::quiet_NaN());
1590 std::vector<double> spdcol(ntrial, std::numeric_limits<double>::quiet_NaN());
1591
1592 for (size_t i = 0; i < ntrial; ++i) {
1593 // dir as string
1594 if (dir_idx < static_cast<int>(ncond)) {
1595 matvar_t *c = cell_at(i, static_cast<size_t>(dir_idx));
1596 if (c && c->class_type == MAT_C_CHAR && c->data) {
1597 size_t len = c->nbytes / c->data_size;
1598 std::string s(static_cast<const char *>(c->data), len);
1599 while (!s.empty() && (s.back() == '\0' || s.back() == ' ')) s.pop_back();
1600 dircol[i] = s;
1601 }
1602 }
1603 // nf numeric
1604 if (nf_idx < static_cast<int>(ncond)) {
1605 matvar_t *c = cell_at(i, static_cast<size_t>(nf_idx));
1606 if (c && c->data) {
1607 if (c->class_type == MAT_C_DOUBLE) {
1608 nfcol[i] = static_cast<const double *>(c->data)[0];
1609 } else if (c->class_type == MAT_C_SINGLE) {
1610 nfcol[i] = static_cast<const float *>(c->data)[0];
1611 } else if (c->class_type == MAT_C_INT32) {
1612 nfcol[i] = static_cast<const int32_t *>(c->data)[0];
1613 }
1614 }
1615 }
1616 // spddxl numeric
1617 if (spd_idx < static_cast<int>(ncond)) {
1618 matvar_t *c = cell_at(i, static_cast<size_t>(spd_idx));
1619 if (c && c->data) {
1620 if (c->class_type == MAT_C_DOUBLE) {
1621 spdcol[i] = static_cast<const double *>(c->data)[0];
1622 } else if (c->class_type == MAT_C_SINGLE) {
1623 spdcol[i] = static_cast<const float *>(c->data)[0];
1624 } else if (c->class_type == MAT_C_INT32) {
1625 spdcol[i] = static_cast<const int32_t *>(c->data)[0];
1626 }
1627 }
1628 }
1629 }
1630
1631 // Subject numeric id
1632 int subj_num = 0;
1633 {
1634 // extract digits from subject_id
1635 for (char ch : cfg.subject_id) {
1636 if (std::isdigit(static_cast<unsigned char>(ch))) {
1637 subj_num = subj_num * 10 + (ch - '0');
1638 }
1639 }
1640 }
1641
1642 // Correction if subject < 8
1643 if (subj_num < 8 && ntrial > 1) {
1644 size_t half = ntrial / 2;
1645 for (size_t i = 0; i < ntrial; ++i) {
1646 spdcol[i] = (i < half) ? 0.04 : 0.08;
1647 }
1648 }
1649
1650 // Build trial indices 1..Ntrial (Matlab-style) based on filters
1651
1652 bool is_loading = (cfg.phase_id == "loading");
1653 std::vector<int> trialnum;
1654 trialnum.reserve(ntrial);
1655 for (size_t i = 0; i < ntrial; ++i) trialnum.push_back(static_cast<int>(i + 1));
1656
1657 for (size_t ii = 0; ii < cfg.nfcond_set.size(); ++ii) {
1658 int nf_set = cfg.nfcond_set[ii];
1659 for (size_t jj = 0; jj < cfg.spddxlcond_set.size(); ++jj) {
1660 double spd_set = cfg.spddxlcond_set[jj];
1661
1662 for (size_t i = 0; i < ntrial; ++i) {
1663 bool pass = true;
1664 if (is_loading) {
1665 bool dir_ok = (dircol[i] == "Ubnf" || dircol[i] == "Rbnf");
1666 bool nf_ok = std::fabs(nfcol[i] - nf_set) < 1e-6;
1667 bool spd_ok = std::fabs(spdcol[i] - spd_set) < 1e-9;
1668 pass = dir_ok && nf_ok && spd_ok;
1669 }
1670 if (pass) trials.push_back(trialnum[i]);
1671 }
1672 }
1673 }
1674
1675 // Clean up
1676 if (cond) Mat_VarFree(cond);
1677 Mat_Close(matfp);
1678
1679 return trials;
1680}
1681
1682bool DicAnalysis::dic2DAnalysis(const std::vector<int>& trial_target, const StagePlan& plan) {
1683 LOG_INFO << "Starting 2D DIC Analysis (using StepDWorkflow)...";
1684
1685 try {
1686 // Create workflow instance
1687 StepDWorkflow workflow(config_);
1688
1689 // Process each trial and stereo pair
1690 for (int trial : trial_target) {
1691 for (int pair = 1; pair <= config_.num_pair; ++pair) {
1692 // Optional narrowing for SLURM array tasks: a single stereopair.
1693 if (plan.only_pair != 0 && pair != plan.only_pair) {
1694 continue;
1695 }
1696 // Format trial as 3-digit string (e.g., "005")
1697 std::ostringstream trial_str;
1698 trial_str << std::setw(3) << std::setfill('0') << trial;
1699
1700 LOG_INFO << "\n========================================";
1701 LOG_INFO << "Processing Trial " << trial << ", Pair " << pair;
1702 LOG_INFO << "========================================";
1703
1704 // Execute workflow (only the sub-steps selected by plan)
1705 auto [outputPath, pairOrder, pairForced] = workflow.execute(trial_str.str(), pair, plan);
1706
1707 if (outputPath.empty()) {
1708 LOG_ERROR << "Workflow failed for trial " << trial << ", pair " << pair;
1709 return false;
1710 }
1711
1712 LOG_INFO << "✓ Trial " << trial << ", pair " << pair << " completed successfully.";
1713 LOG_INFO << " Output path: " << outputPath;
1714 LOG_INFO << " Pair order: [" << pairOrder[0] << ", " << pairOrder[1] << "]";
1715 LOG_INFO << " Pair forced: " << (pairForced ? "true" : "false");
1716 }
1717 }
1718
1719 LOG_INFO << "\n✓ All 2D DIC Analysis completed successfully!";
1720 return true;
1721
1722 } catch (const std::exception& e) {
1723 LOG_ERROR << "in 2D DIC Analysis: " << e.what();
1724 return false;
1725 }
1726}
1727
1728
1729std::vector<std::string> DicAnalysis::loadImageSequence(const std::string& trial_path) {
1730 std::vector<std::string> images;
1731
1732 // Look for common image formats
1733 std::vector<std::string> extensions = {".png", ".jpg", ".jpeg", ".tiff", ".bmp"};
1734
1735 try {
1736 if (std::filesystem::exists(trial_path)) {
1737 for (const auto& entry : std::filesystem::directory_iterator(trial_path)) {
1738 if (entry.is_regular_file()) {
1739 std::string filename = entry.path().string();
1740 for (const auto& ext : extensions) {
1741 if (filename.size() >= ext.size() &&
1742 filename.compare(filename.size() - ext.size(), ext.size(), ext) == 0) {
1743 images.push_back(filename);
1744 break;
1745 }
1746 }
1747 }
1748 }
1749 }
1750
1751 // Sort images to ensure proper sequence
1752 std::sort(images.begin(), images.end());
1753
1754 } catch (const std::exception& e) {
1755 LOG_ERROR << "loading image sequence: " << e.what();
1756 }
1757
1758 return images;
1759}
1760
1761bool DicAnalysis::setupNcorrAnalysis(const std::vector<std::string>& images,
1762 DIC_analysis_input& dic_input) {
1763 try {
1764 if (images.size() < 2) {
1765 LOG_ERROR << "Need at least 2 images for DIC analysis";
1766 return false;
1767 }
1768
1769 // Convert string paths to Image2D objects
1770 std::vector<Image2D> ncorr_images;
1771 for (const auto& img_path : images) {
1772 ncorr_images.emplace_back(img_path);
1773 }
1774
1775 Image2D first_image(images[0]);
1776 ROI2D roi(first_image.get_gs() > 0.1);
1777
1778 // Setup DIC input with parameters similar to the ncorr test
1779 dic_input = DIC_analysis_input(
1780 ncorr_images, // Images
1781 roi, // ROI
1782 3, // scalefactor
1783 INTERP::QUINTIC_BSPLINE_PRECOMPUTE, // Interpolation
1784 SUBREGION::CIRCLE, // Subregion shape
1785 config_.subregion_radius, // Subregion radius
1786 4, // # of threads
1787 DIC_analysis_config::NO_UPDATE, // DIC configuration
1788 config_.debug_mode // Debugging enabled/disabled
1789 );
1790
1791 return true;
1792
1793 } catch (const std::exception& e) {
1794 LOG_ERROR << "setting up ncorr analysis: " << e.what();
1795 return false;
1796 }
1797}
int ref_trial_id
Definition config.h:75
bool generate_videos
Definition config.h:228
bool step_e_distortion_removal
Definition config.h:157
void overrideSubject(const std::string &subject)
Definition config.cpp:622
bool mapLogic
Definition config.h:91
bool step_f_compute_rbm
Definition config.h:165
double vid_sample_freq
Definition config.h:43
bool step_e_ray_bvh
Definition config.h:156
std::string subject_id
Definition config.h:65
std::string step_e_stitch_mode
Definition config.h:151
void updateVariables()
Definition config.cpp:18
double step_f_freq_filt
Definition config.h:164
std::vector< std::pair< int, int > > camera_pairs
Definition config.h:37
bool step_f_temporal_filtering
Definition config.h:163
int num_pair
Definition config.h:34
std::vector< std::string > plotopt
Definition config.h:92
bool debug_mode
Definition config.h:84
int subregion_radius
Definition config.h:131
std::string phase_id
Definition config.h:66
std::string fileversion
Definition config.h:125
std::vector< int > nfcond_set
Definition config.h:68
std::vector< double > spddxlcond_set
Definition config.h:69
std::string data_format
Definition config.h:168
DicAnalysis(const Config &config)
std::vector< int > searchTrialTarget(const std::string &subject)
Search the protocol .mat for the list of trial IDs of a given subject.
bool run()
Run complete DIC analysis pipeline Executes preprocessing, Step D, Step E, and Step F in sequence.
static void setCameraPairs(const std::vector< std::pair< int, int > > &camera_pairs)
Register the process-wide stereopair -> (cam_first, cam_second) table.
Definition utils.cpp:246
static std::vector< std::string > findFiles(const std::string &directory, const std::string &pattern)
Definition utils.cpp:193
static std::string buildDic3DPPresultsFilePath(const std::string pathdir, int num_pair, const std::string &fileversion, const std::string &ext)
Definition utils.cpp:1346
static std::string buildDic2DPairResultsFilePath(const std::string pathdir, int cam_1, int cam_2, const std::string &extension)
Definition utils.cpp:1357
static std::string buildProtocolDir(const Config &config, bool with_data_or_dic_path=false, bool with_rawdata=false, bool with_speckles=false, bool with_material=false)
Build protocol path.
Definition utils.cpp:1241
static std::string buildCalibDir(const Config &config)
Build the calib path.
Definition utils.cpp:1257
static std::string buildVizPath(const std::string pathdir, const std::string filename)
Definition utils.cpp:1353
static void undistortPoints(const std::vector< cv::Point2d > &points_in, const CameraParameters &params, std::vector< cv::Point2d > &points_out)
Undistort 2D points using camera parameters Matches MATLAB: undistortPoints(points,...
Definition utils.cpp:911
static bool computeRigidTransform(const std::vector< Eigen::Vector3d > &points_from, const std::vector< Eigen::Vector3d > &points_to, RigidTransform &transform)
Compute optimal rigid body transformation between two point clouds Matches MATLAB: rigidTransformatio...
Definition utils.cpp:1060
static bool DLT11Calibration(const double *P2, const double *P3, size_t N, std::vector< double > &L)
Compute DLT (Direct Linear Transformation) 11 parameters from 2D-3D point correspondences Matches MAT...
Definition utils.cpp:999
static std::string buildOutputUntilPhaseDir(const Config &config, int trial)
Build output dir until the Phase folder.
Definition utils.cpp:1195
static std::vector< Eigen::Vector3d > applyRigidTransform(const std::vector< Eigen::Vector3d > &points_in, const RigidTransform &transform)
Apply rigid transformation to point cloud points_out = R * points_in + t.
Definition utils.cpp:1140
static bool loadCameraParameters(const std::string &mat_path, CameraParameters &params)
Load camera calibration parameters from MATLAB .mat file Matches MATLAB: load('cameraCBparameters_cam...
Definition utils.cpp:804
static std::string buildDic3DCombinedFilePath(const std::string pathdir, int num_pair, const std::string &ext)
Definition utils.cpp:1341
static void getCamerasForPair(int stereopair, int &cam_first, int &cam_second)
Get camera numbers for a given stereopair using the registered table (setCameraPairs).
Definition utils.cpp:254
static std::unique_ptr< DataSerializer > create(const std::string &format)
Create a serializer for the given format string.
static bool loadDLTCalibration(const std::string &mat_path, DLTCalibrationData &calib)
static bool loadProtocol(const std::string &protocol_path, ProtocolFileData &protocol)
StepDWorkflow class Implements the complete 2D DIC analysis workflow Equivalent to stepD_2DDIC....
Visualization manager class Handles loading, filtering, and exporting DIC results.
Lightweight, dependency-free logging facility for CPPxDIC.
#define LOG_DEBUG
Definition logging.h:169
#define LOG_ERROR
Definition logging.h:172
#define LOG_WARN
Definition logging.h:171
#define LOG_INFO
Definition logging.h:170
DataSerializer - Abstract base class for format-agnostic DIC data I/O.
DIC3Dcombined stitchPairsSimple(const std::vector< DIC3DpairResults > &all_pairs)
Simple append-based stitching for multiple stereo pairs.
std::tuple< std::vector< std::vector< double > >, std::vector< std::vector< double > >, std::vector< std::vector< double > > > filterTime3D(const std::vector< std::vector< double > > &disp_x, const std::vector< std::vector< double > > &disp_y, const std::vector< std::vector< double > > &disp_z, double freq_filt=10.0, double freq_acq=50.0)
Apply temporal filtering to 3D displacement vectors.
void setStitchRayAccel(bool enabled)
Enable/disable the AABB-hierarchy prefilter used by the overlap-removal ray casts.
std::vector< double > computeFaceIsotropyIndex(const std::vector< int > &faces, const std::vector< Eigen::Vector3d > &vertices)
Compute isotropy index for triangular faces.
FrameDeformationResult computeTriSurfaceDeformation(const std::vector< int > &faces, const std::vector< Eigen::Vector3d > &vertices_ref, const std::vector< std::vector< Eigen::Vector3d > > &vertices_def, bool cumulative=true)
Compute deformation and strain for triangular surface elements.
DIC3Dcombined stitchPairsGeometric(const std::vector< DIC3DpairResults > &all_pairs, const std::vector< int > &pair_order, bool pair_forced=false)
Geometric stitching with overlap removal and boundary zipping Matches MATLAB's DIC3DsurfaceStitch....
StagePlan — selects which pipeline sub-steps a single process executes.
Definition stage_plan.h:21
int only_pair
0 = all stereopairs; otherwise restrict to this pair
Definition stage_plan.h:31
bool match
ROI/seed + REF->trial matching + camera-to-camera matching.
Definition stage_plan.h:23
bool deform
Step F: deformation / strain.
Definition stage_plan.h:28
bool recon
Step E: 3D reconstruction.
Definition stage_plan.h:27
bool format
combine ncorr outputs -> myDIC2DpairResults
Definition stage_plan.h:25
bool anyStepD() const
True if any Step-D sub-step is requested.
Definition stage_plan.h:35
bool track
per-camera NCorr tracking
Definition stage_plan.h:24
Eigen::Vector3d t
Definition utils.h:114
Eigen::Matrix3d R
Definition utils.h:113
DIC 2D Pair Results Equivalent to MATLAB's DIC2DpairResults.
std::vector< Points2D > Points
std::vector< std::vector< double > > CorCoeffVec
std::vector< double > FaceColors
std::vector< int > Faces
DIC 3D Post-Processing Results Equivalent to MATLAB's DIC3DPPresults.
std::vector< std::vector< double > > FaceIsoInd
std::vector< std::vector< double > > Points3D_ARBM_y
std::vector< std::vector< double > > Points3D_ARBM_x
FrameDeformationResult deform_full
std::vector< std::vector< double > > Points3D_ARBM_z
DIC 3D Combined Results Equivalent to MATLAB's DIC3Dcombined This represents stitched results from mu...
std::vector< int > FacePairInds
DistortionData distortion
std::vector< int > PointPairInds
std::vector< DIC3DpairResults > AllPairsResults
std::vector< Point3DFrame > Points3D
std::vector< int > pairIndices
std::vector< double > FaceColors
std::vector< std::vector< double > > FaceCorrComb
std::vector< std::vector< double > > corrComb
std::vector< DIC2DPairResults > DIC2Dinfo
CalibrationData calibration
std::vector< int > Faces
std::vector< std::vector< double > > FaceCentroids
DIC 3D Pair Results (Individual stereo pair) Equivalent to MATLAB's DIC3DpairResults This represents ...
std::vector< std::vector< double > > DLTparameters
std::vector< std::vector< double > > corrComb
std::vector< int > Faces
std::vector< std::vector< double > > FaceCentroids
std::vector< double > FaceColors
std::vector< Point3DFrame > Points3D
std::vector< std::string > DLTpath
std::vector< std::vector< double > > FaceCorrComb
std::vector< std::string > distortionPath
std::vector< std::string > distortionModel
std::vector< int > cameraPairInd
std::vector< DeformGradient > F
std::vector< StrainTensor > strain
std::vector< std::vector< double > > princStrain
std::vector< std::vector< double > > maxShearStrain
Deformation gradient tensor.
std::vector< double > F13
std::vector< double > F23
std::vector< double > F12
std::vector< double > F22
std::vector< double > F21
std::vector< double > F32
std::vector< double > F31
std::vector< double > F33
std::vector< double > F11
std::vector< std::vector< double > > DispVec
std::vector< std::vector< double > > DispMgn
Frame-based deformation result Contains deformation results for all frames.
std::vector< DeformationResult > frames
3D Points per frame
std::vector< double > z
std::vector< double > x
std::vector< double > y
std::vector< TrialEntry > trials
Definition mat_reader.h:43
std::vector< std::vector< double > > TransVec
std::vector< std::vector< double > > RotMat
std::vector< double > E11
std::vector< double > E12
std::vector< double > E13
std::vector< double > E21
std::vector< double > E31
std::vector< double > E32
std::vector< double > E23
std::vector< double > E22
std::vector< double > E33