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";
45 for (
int trial : trial_target) {
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!";
59 LOG_INFO <<
"Loading DIC3Dcombined from: " << dic3d_file;
62 if (!serializer->loadDIC3Dcombined(dic3d_file, dic3d)) {
63 LOG_ERROR <<
"Failed to load DIC3Dcombined structure";
68 LOG_ERROR <<
"DIC3Dcombined has no 3D data";
73 << dic3d.
Points3D[0].x.size() <<
" points, "
74 << dic3d.
Faces.size() / 3 <<
" faces";
76 << (dic3d.
corrComb.empty() ?
"" :
" (" + std::to_string(dic3d.
corrComb[0].size()) +
" per frame)");
84 size_t inf_pts_total = 0;
85 size_t bad_face_total = 0;
87 for (
size_t f = 0; f < std::min(dic3d.
Points3D.size(),
size_t(3)); ++f) {
89 for (
size_t i = 0; i < nP; ++i) {
95 LOG_DEBUG <<
" [diag] Frame " << f <<
": " << nan_pts <<
"/" << nP <<
" NaN points";
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)
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!";
110 LOG_INFO <<
"\nConverting data to Eigen format...";
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(
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(
138 vertices_all_frames.push_back(std::move(frame_verts));
141 LOG_INFO <<
" Converted " << vertices_all_frames.size() <<
" frames";
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) {
153 counts.push_back(nan_count);
158 auto count_points_with_any_nan = [](
const std::vector<std::vector<Eigen::Vector3d>>& frames) {
159 if (frames.empty()) {
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()) {
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) {
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()) {
194 if (verts[v0].hasNaN() || verts[v1].hasNaN() || verts[v2].hasNaN()) {
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()) {
211 max_disp = std::max(max_disp, (frame[ipt] - ref[ipt]).norm());
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
222 if (!raw_nan_counts.empty()) {
223 LOG_DEBUG <<
" Raw NaN points in frame 0: " << raw_nan_counts[0]
226 LOG_DEBUG <<
" Raw valid faces in frame 0: "
227 << count_valid_faces(vertices_all_frames[0])
229 LOG_DEBUG <<
" Raw max displacement magnitude: " << raw_max_disp;
232 bool used_temporal_filtering =
false;
234 LOG_INFO <<
"\nApplying temporal filtering...";
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));
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();
253 auto [filt_x, filt_y, filt_z] =
filterTime3D(disp_x, disp_y, disp_z, freq_filt, freq_acq);
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];
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;
273 LOG_DEBUG <<
" Filtered points with any NaN over time: " << filtered_any_nan_points
275 if (!filtered_nan_counts.empty()) {
276 LOG_DEBUG <<
" Filtered NaN points in frame 0: " << filtered_nan_counts[0]
279 LOG_DEBUG <<
" Filtered valid faces in frame 0: "
280 << filtered_valid_faces <<
"/" << nFaces;
281 LOG_DEBUG <<
" Filtered max displacement magnitude: " << filtered_max_disp;
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";
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];
296 used_temporal_filtering =
true;
297 LOG_INFO <<
" ✓ Temporal filtering applied (freqFilt=" << freq_filt
298 <<
" Hz, freqAcq=" << freq_acq <<
" Hz)";
301 LOG_INFO <<
"\nSkipping temporal filtering (too few frames: "
302 << vertices_all_frames.size() <<
")";
304 LOG_INFO <<
"\nTemporal filtering disabled (step_f_temporal_filtering=false)";
308 LOG_INFO <<
"\nUpdating Points3D with filtered data...";
310 dic3d.
Points3D.resize(vertices_all_frames.size());
311 for (
size_t iframe = 0; iframe < vertices_all_frames.size(); ++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();
322 if (used_temporal_filtering) {
323 LOG_INFO <<
" ✓ Points3D updated with filtered data";
325 LOG_INFO <<
" ✓ Points3D kept from unfiltered geometry";
329 LOG_INFO <<
"\nRecomputing displacement after filtering...";
332 dic3d.
Disp.
DispVec.resize(vertices_all_frames.size());
333 dic3d.
Disp.
DispMgn.resize(vertices_all_frames.size());
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];
339 disp_vec.resize(nPoints * 3);
340 disp_mgn.resize(nPoints);
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();
347 disp_vec[ipt * 3 + 0] = dx;
348 disp_vec[ipt * 3 + 1] = dy;
349 disp_vec[ipt * 3 + 2] = dz;
351 disp_mgn[ipt] = std::sqrt(dx*dx + dy*dy + dz*dz);
354 LOG_INFO <<
" ✓ Displacement recomputed for " << vertices_all_frames.size() <<
" frames";
357 LOG_INFO <<
"\nRecomputing face centroids after filtering...";
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);
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];
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;
384 centroids[iface * 3 + 0] = cx;
385 centroids[iface * 3 + 1] = cy;
386 centroids[iface * 3 + 2] = cz;
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();
394 LOG_INFO <<
" ✓ Face centroids recomputed for " << vertices_all_frames.size() <<
" frames";
398 LOG_INFO <<
"\nRecomputing face correlation after filtering...";
399 if (!dic3d.
corrComb.empty() && dic3d.
corrComb.size() >= vertices_all_frames.size()) {
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);
407 const auto& point_corr = dic3d.
corrComb[iframe];
408 size_t nCorr = point_corr.size();
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];
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]});
420 face_corr[iface] = std::numeric_limits<double>::quiet_NaN();
424 LOG_INFO <<
" ✓ Face correlation recomputed for " << vertices_all_frames.size() <<
" frames";
426 LOG_WARN <<
" ⚠ corrComb not available (empty or frame count mismatch), skipping face correlation";
430 std::vector<Utils::RigidTransform> rbm_transforms;
431 std::vector<std::vector<Eigen::Vector3d>> vertices_all_frames_ARBM;
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());
438 for (
size_t iframe = 0; iframe < vertices_all_frames.size(); ++iframe) {
441 vertices_all_frames[iframe],
446 rbm_transforms[iframe] = transform;
450 vertices_all_frames_ARBM.push_back(verts_arbm);
452 LOG_WARN <<
" RBM computation failed for frame " << iframe;
453 vertices_all_frames_ARBM.push_back(vertices_all_frames[iframe]);
456 LOG_INFO <<
" ✓ RBM transformations computed for " << vertices_all_frames.size() <<
" frames";
458 LOG_INFO <<
"\nRBM/ARBM computation disabled (step_f_compute_rbm=false)";
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)";
473 LOG_INFO <<
" ✓ Deformation computation complete";
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)
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";
497 LOG_INFO <<
"\nComputing 3D surface deformation (after RBM removal)...";
500 vertices_all_frames_ARBM[0],
501 vertices_all_frames_ARBM,
504 LOG_INFO <<
" ✓ Deformation computation complete (ARBM)";
508 LOG_INFO <<
"\nBuilding DIC3DPPresults structure...";
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) {
520 vertices_all_frames[iframe]
523 LOG_INFO <<
" ✓ Face isotropy index computed for " << vertices_all_frames.size() <<
" frames";
554 for (
size_t iframe = 0; iframe < deform_result.
n_frames; ++iframe) {
555 const auto& frame = deform_result.
frames[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);
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);
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);
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);
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];
615 maxShear = frame.EShearMax;
618 LOG_INFO <<
" Populated deformation data for " << ppresults.
n_frames <<
" frames";
622 LOG_INFO <<
" Populating ARBM deformation data...";
628 for (
size_t iframe = 0; iframe < deform_result_ARBM.
n_frames; ++iframe) {
629 const auto& frame = deform_result_ARBM.
frames[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);
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);
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);
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);
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];
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];
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)};
690 for (
size_t iframe = 0; iframe < vertices_all_frames_ARBM.size(); ++iframe) {
691 const auto& verts = vertices_all_frames_ARBM[iframe];
695 for (
size_t i = 0; i < verts.size(); ++i) {
702 LOG_INFO <<
" ✓ Populated all deformation data (with RBM and ARBM)";
710 if (std::filesystem::exists(ppout)) {
711 LOG_INFO <<
"\nCheckpoint found: " << ppout <<
" (skipping)";
713 LOG_INFO <<
"\nWriting DIC3DPPresults to: " << ppout;
714 if (!pp_serializer->saveDIC3DPPresults(ppout, ppresults)) {
715 LOG_ERROR <<
"Failed to write DIC3DPPresults";
718 LOG_INFO <<
"✓ Saved DIC3DPPresults: " << ppout;
724 LOG_INFO <<
"\n=== Generating Visualization Exports ===";
729 viz.printTrialInfo(ppresults);
732 viz.applyTemporalFilter(ppresults);
735 auto vis_data = viz.prepareVisualizationData(ppresults);
738 std::filesystem::create_directories(output_dir +
"viz");
742 viz.exportData(vis_data, viz_path);
746 viz.generateSummaryStats(ppresults, stats_path);
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);
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");
764 LOG_INFO <<
"✓ Visualization exports complete";
766 }
catch (
const std::exception& e) {
767 LOG_WARN <<
"Visualization export failed: " << e.what();
768 LOG_WARN <<
" (Continuing anyway...)";
772 LOG_INFO <<
"\n=== Step F Complete ===";
773 LOG_INFO <<
"✓ 3D deformation and strain analysis finished for trial " << trial;
778 }
catch (
const std::exception& e) {
779 LOG_ERROR <<
"in Deformation/Strain Analysis: " << e.what();
783bool DicAnalysis::dic3DReconstruction(
const std::vector<int>& trial_target) {
784 LOG_INFO <<
"Starting 3D Reconstruction (Step E)...";
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;
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];
805 Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 4, 3>> qr(M);
807 double nan = std::numeric_limits<double>::quiet_NaN();
808 return std::array<double, 3>{nan, nan, nan};
811 Eigen::Matrix<double, 3, 1> X = qr.solve(V);
812 return std::array<double, 3>{X[0], X[1], X[2]};
815 for (
int trial : trial_target) {
817 std::vector<DIC3DpairResults> all_pairs;
818 std::vector<DIC2DPairResults> dic2d_info;
821 std::ostringstream trial_str;
822 trial_str << std::setw(3) << std::setfill(
'0') << trial;
831 std::set<int> unique_cams;
832 for (
int p = 1; p <= config_.
num_pair; ++p) {
835 unique_cams.insert(c1);
836 unique_cams.insert(c2);
840 std::map<int, DLTCalibrationData> dlt_all_cams;
841 for (
int cam_id : unique_cams) {
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();
857 if (dlt_path.empty()) {
858 LOG_ERROR <<
"DLT calibration file not found for camera " << cam_id
859 <<
" in " << calib_dir;
865 dlt_all_cams[cam_id] = std::move(calib);
867 LOG_ERROR <<
"Failed to load DLT calibration for camera " << cam_id;
871 if (dlt_all_cams.empty()) {
872 LOG_ERROR <<
"No DLT calibrations loaded. Cannot proceed with 3D reconstruction.";
875 LOG_INFO <<
"Loaded DLT calibrations for " << dlt_all_cams.size() <<
" cameras";
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) {
894 distortion_all_cams[cam_id] = params;
895 distortion_paths[cam_id] = entry.path().string();
903 if (use_distortion_removal) {
904 LOG_INFO <<
"Distortion removal enabled for all " << unique_cams.size() <<
" cameras";
915 for (
int cam_id : unique_cams) {
916 auto& dlt = dlt_all_cams[cam_id];
917 const auto& dist = distortion_all_cams[cam_id];
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";
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]);
932 std::vector<cv::Point2d> centroids_undist;
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;
945 size_t d0 = dlt.C3Dtrue_dim0;
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;
950 for (
size_t row = 0; row < d0; ++row) {
951 size_t dst_idx = (ci * d0 + row);
952 size_t src_idx = (row * dlt.C3Dtrue_dim1 + col);
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];
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";
967 std::vector<double> L_new;
969 dlt.DLTparams = L_new;
970 LOG_INFO <<
"Camera " << cam_id <<
": DLT params recalculated from "
971 << nPoints <<
" undistorted calibration points";
973 LOG_ERROR <<
"Camera " << cam_id <<
": DLT11Calibration failed";
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;
985 for (
int pair = 1; pair <= config_.
num_pair; ++pair) {
986 LOG_INFO <<
"\n=== Processing Pair " << pair <<
" ===";
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.";
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;
1013 if (!std::filesystem::exists(dic2d_file)) {
1014 LOG_ERROR <<
"Missing DIC2DPairResults for trial " << trial <<
", pair " << pair <<
". Skipping.";
1015 LOG_ERROR <<
" Expected: " << dic2d_file;
1020 if (!d_serializer->loadDIC2DPairResults(dic2d_file, DIC2D)) {
1021 LOG_ERROR <<
"Failed to load DIC2DPairResults: " << dic2d_file <<
". Skipping.";
1026 LOG_DEBUG <<
"Loaded DIC2DPairResults from: " << dic2d_file;
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;
1036 const auto& F = DIC2D.
Faces;
1038 const auto& Points = DIC2D.
Points;
1043 pair_result.
DLTpath = {dlt_cam1.filePath, dlt_cam2.filePath};
1045 pair_result.
Faces = F;
1049 if (use_distortion_removal) {
1051 pair_result.
distortionPath = {distortion_paths[cam_1], distortion_paths[cam_2]};
1057 size_t nFaces = F.size() / 3;
1058 std::vector<double> P3D_ref;
1061 for (
int ii = 0; ii < nImages; ++ii) {
1063 const auto& P1 = Points[ii];
1064 const auto& P2 = Points[ii + nImages];
1065 size_t N = P1.x.size();
1068 std::vector<double> pts3d(N * 3, std::numeric_limits<double>::quiet_NaN());
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];
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();
1084 pair_result.
Points3D.push_back(frame_pts);
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);
1095 pair_result.
corrComb.push_back(corr_comb);
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]});
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;
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);
1135 all_pairs.push_back(pair_result);
1136 dic2d_info.push_back(std::move(DIC2D));
1138 LOG_INFO <<
"✓ Pair " << pair <<
" complete: "
1139 << pair_result.
Points3D[0].x.size() <<
" points, "
1140 << nFaces <<
" faces";
1144 LOG_INFO <<
"\n=== Stitching " << all_pairs.size() <<
" pairs ===";
1146 if (all_pairs.empty()) {
1147 LOG_ERROR <<
"No pairs successfully reconstructed for trial " << trial;
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;
1155 for (
const auto& dic2d : dic2d_info) {
1156 if (dic2d.pairOrder.empty()) {
1157 LOG_ERROR <<
"Empty pairOrder in DIC2D for trial " << trial;
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;
1174 if (stitch_metadata_mismatch) {
1175 LOG_WARN <<
"Skipping trial " << trial <<
" due to stitch metadata mismatch";
1181 LOG_INFO <<
" step_e_stitch_mode=simple: appending " << all_pairs.size()
1182 <<
" pairs without overlap removal";
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;
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);
1197 if (!protocol_files.empty()) {
1198 protocol_loaded =
true;
1202 protocol_available =
false;
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;
1215 stitch_pair_forced =
false;
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);
1223 LOG_INFO <<
" Stitch metadata missing from DIC2D results; using legacy fallback";
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];
1232 LOG_INFO <<
" Stitch order: " << stitch_order_oss.str();
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)");
1238 LOG_INFO <<
"Geometric Stitching done!";
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";
1262 if (std::filesystem::exists(e_out)) {
1263 LOG_INFO <<
"Checkpoint found: " << e_out <<
" (skipping)";
1265 LOG_INFO <<
"Writing DIC3Dcombined to: " << e_out;
1270 const std::string e_tmp = e_out +
".partial";
1271 if (!e_serializer->saveDIC3Dcombined(e_tmp, stitched)) {
1272 LOG_ERROR <<
"Failed to write DIC3Dcombined";
1274 std::filesystem::remove(e_tmp, ec);
1277 std::filesystem::rename(e_tmp, e_out, ec);
1279 LOG_ERROR <<
"Failed to finalize DIC3Dcombined: " << ec.message();
1281 LOG_INFO <<
"✓ Saved DIC3Dcombined: " << e_out;
1288 }
catch (
const std::exception& e) {
1289 LOG_ERROR <<
"in 3D Reconstruction: " << e.what();
1294bool DicAnalysis::setupNcorrAnalysis(
const std::vector<std::string>& images,
1295 const std::string& roi_mask_image_path,
1296 DIC_analysis_input& dic_input) {
1298 if (images.size() < 2) {
1299 LOG_ERROR <<
"Need at least 2 images for DIC analysis";
1303 std::vector<Image2D> ncorr_images;
1304 for (
const auto& img_path : images) {
1305 ncorr_images.emplace_back(img_path);
1309 Image2D roi_mask(roi_mask_image_path);
1310 ROI2D roi(roi_mask.get_gs() > 0.5);
1312 dic_input = DIC_analysis_input(
1316 INTERP::QUINTIC_BSPLINE_PRECOMPUTE,
1320 DIC_analysis_config::NO_UPDATE,
1326 }
catch (
const std::exception& e) {
1327 LOG_ERROR <<
"setting up ncorr analysis with ROI mask: " << e.what();
1335 std::vector<int> trial_target = {7};
1336 return run(trial_target);
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 <<
", ";
1352 trial_target_oss <<
"]";
1353 LOG_INFO <<
"Trial target set: " << trial_target_oss.str();
1358 std::string ext = serializer->extension();
1361 auto check_2d_outputs_exist = [&]() ->
bool {
1362 for (
int trial : trial_target) {
1364 for (
int pair = 1; pair <= config_.
num_pair; ++pair) {
1369 if (!std::filesystem::exists(path)) {
1378 auto check_3d_outputs_exist = [&]() ->
bool {
1379 for (
int trial : trial_target) {
1382 if (!std::filesystem::exists(dic3d_file)) {
1390 auto check_deformation_outputs_exist = [&]() ->
bool {
1391 for (
int trial : trial_target) {
1394 if (!std::filesystem::exists(pp_file)) {
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)";
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);
1423 LOG_INFO <<
"✓ DIC 2D Analysis done in " << duration.count() / 1000.0 <<
" s";
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)";
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);
1443 LOG_ERROR <<
"3D Reconstruction failed!";
1447 LOG_INFO <<
"✓ DIC 3D Reconstruction done in " << duration.count() / 1000.0 <<
" s";
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)";
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);
1467 LOG_ERROR <<
"Deformation Analysis failed!";
1471 LOG_INFO <<
"✓ DIC Deformation Analysis done in " << duration.count() / 1000.0 <<
" s";
1475 LOG_INFO <<
"\n========================================";
1476 LOG_INFO <<
"✓ All DIC analysis steps complete!";
1477 LOG_INFO <<
"========================================";
1482std::vector<int> DicAnalysis::searchTrialTarget() {
1484 return searchTrialTarget(config_.
subject_id);
1488 std::vector<int> trials;
1501 if (protos.empty()) {
1502 LOG_ERROR <<
"Protocol file not found in: " << protocol_dir;
1509 std::string proto_file = protos.front();
1512 mat_t *matfp = Mat_Open(proto_file.c_str(), MAT_ACC_RDONLY);
1514 LOG_ERROR <<
"Failed to open MAT file: " << proto_file;
1521 matvar_t *cond = Mat_VarRead(matfp,
"cond");
1522 if (!cond || cond->class_type != MAT_C_STRUCT) {
1523 if (cond) Mat_VarFree(cond);
1525 LOG_ERROR <<
"Variable 'cond' not found or not a struct in: " << proto_file;
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);
1537 LOG_ERROR <<
"Fields 'titles' or 'table' missing or of wrong type in 'cond'";
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];
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);
1554 while (!name.empty() && (name.back() ==
'\0' || name.back() ==
' ')) name.pop_back();
1556 col_names.push_back(name);
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);
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);
1572 LOG_ERROR <<
"Required columns not found in titles (need 'dir','nf','spddxl')";
1579 size_t ntrial = table->dims[0];
1580 size_t ncond = table->dims[1];
1583 auto cell_at = [&](
size_t i,
size_t j) -> matvar_t * {
1584 size_t idx = i + j * ntrial;
1585 return static_cast<matvar_t **
>(table->data)[idx];
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());
1592 for (
size_t i = 0; i < ntrial; ++i) {
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();
1604 if (nf_idx <
static_cast<int>(ncond)) {
1605 matvar_t *c = cell_at(i,
static_cast<size_t>(nf_idx));
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];
1617 if (spd_idx <
static_cast<int>(ncond)) {
1618 matvar_t *c = cell_at(i,
static_cast<size_t>(spd_idx));
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];
1636 if (std::isdigit(
static_cast<unsigned char>(ch))) {
1637 subj_num = subj_num * 10 + (ch -
'0');
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;
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));
1657 for (
size_t ii = 0; ii < cfg.
nfcond_set.size(); ++ii) {
1662 for (
size_t i = 0; i < ntrial; ++i) {
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;
1670 if (pass) trials.push_back(trialnum[i]);
1676 if (cond) Mat_VarFree(cond);
1682bool DicAnalysis::dic2DAnalysis(
const std::vector<int>& trial_target,
const StagePlan& plan) {
1683 LOG_INFO <<
"Starting 2D DIC Analysis (using StepDWorkflow)...";
1690 for (
int trial : trial_target) {
1691 for (
int pair = 1; pair <= config_.
num_pair; ++pair) {
1697 std::ostringstream trial_str;
1698 trial_str << std::setw(3) << std::setfill(
'0') << trial;
1700 LOG_INFO <<
"\n========================================";
1701 LOG_INFO <<
"Processing Trial " << trial <<
", Pair " << pair;
1702 LOG_INFO <<
"========================================";
1705 auto [outputPath, pairOrder, pairForced] = workflow.execute(trial_str.str(), pair, plan);
1707 if (outputPath.empty()) {
1708 LOG_ERROR <<
"Workflow failed for trial " << trial <<
", pair " << pair;
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");
1719 LOG_INFO <<
"\n✓ All 2D DIC Analysis completed successfully!";
1722 }
catch (
const std::exception& e) {
1723 LOG_ERROR <<
"in 2D DIC Analysis: " << e.what();
1729std::vector<std::string> DicAnalysis::loadImageSequence(
const std::string& trial_path) {
1730 std::vector<std::string> images;
1733 std::vector<std::string> extensions = {
".png",
".jpg",
".jpeg",
".tiff",
".bmp"};
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);
1752 std::sort(images.begin(), images.end());
1754 }
catch (
const std::exception& e) {
1755 LOG_ERROR <<
"loading image sequence: " << e.what();
1761bool DicAnalysis::setupNcorrAnalysis(
const std::vector<std::string>& images,
1762 DIC_analysis_input& dic_input) {
1764 if (images.size() < 2) {
1765 LOG_ERROR <<
"Need at least 2 images for DIC analysis";
1770 std::vector<Image2D> ncorr_images;
1771 for (
const auto& img_path : images) {
1772 ncorr_images.emplace_back(img_path);
1775 Image2D first_image(images[0]);
1776 ROI2D roi(first_image.get_gs() > 0.1);
1779 dic_input = DIC_analysis_input(
1783 INTERP::QUINTIC_BSPLINE_PRECOMPUTE,
1787 DIC_analysis_config::NO_UPDATE,
1793 }
catch (
const std::exception& e) {
1794 LOG_ERROR <<
"setting up ncorr analysis: " << e.what();
bool step_e_distortion_removal
void overrideSubject(const std::string &subject)
std::string step_e_stitch_mode
std::vector< std::pair< int, int > > camera_pairs
bool step_f_temporal_filtering
std::vector< std::string > plotopt
std::vector< int > nfcond_set
std::vector< double > spddxlcond_set
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.
static std::vector< std::string > findFiles(const std::string &directory, const std::string &pattern)
static std::string buildDic3DPPresultsFilePath(const std::string pathdir, int num_pair, const std::string &fileversion, const std::string &ext)
static std::string buildDic2DPairResultsFilePath(const std::string pathdir, int cam_1, int cam_2, const std::string &extension)
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.
static std::string buildCalibDir(const Config &config)
Build the calib path.
static std::string buildVizPath(const std::string pathdir, const std::string filename)
static void undistortPoints(const std::vector< cv::Point2d > &points_in, const CameraParameters ¶ms, std::vector< cv::Point2d > &points_out)
Undistort 2D points using camera parameters Matches MATLAB: undistortPoints(points,...
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...
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...
static std::string buildOutputUntilPhaseDir(const Config &config, int trial)
Build output dir until the Phase folder.
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.
static bool loadCameraParameters(const std::string &mat_path, CameraParameters ¶ms)
Load camera calibration parameters from MATLAB .mat file Matches MATLAB: load('cameraCBparameters_cam...
static std::string buildDic3DCombinedFilePath(const std::string pathdir, int num_pair, const std::string &ext)
static void getCamerasForPair(int stereopair, int &cam_first, int &cam_second)
Get camera numbers for a given stereopair using the registered table (setCameraPairs).
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.
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.
int only_pair
0 = all stereopairs; otherwise restrict to this pair
bool match
ROI/seed + REF->trial matching + camera-to-camera matching.
bool deform
Step F: deformation / strain.
bool recon
Step E: 3D reconstruction.
bool format
combine ncorr outputs -> myDIC2DpairResults
bool anyStepD() const
True if any Step-D sub-step is requested.
bool track
per-camera NCorr tracking
DIC 2D Pair Results Equivalent to MATLAB's DIC2DpairResults.
std::vector< Points2D > Points
std::vector< std::vector< double > > CorCoeffVec
std::vector< double > FaceColors
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< 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< 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< std::vector< double > > DispVec
std::vector< std::vector< double > > DispMgn
std::vector< TrialEntry > trials
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