fix(lib): correct main.cpp voxel ordering + harden simulateMTR per review

FIX 1: Rebuild spatialCoords in SIMPLNX z,y,x order (iz outer, iy middle,
ix inner) so row k matches sim.phi*/mtrIndex[k]. The old loop used z-x-y
order which mismatched coordinates with orientations for non-square grids.

FIX 2: Replace magic literals 72, 36, 72 in the simulateMTR call with
named constants (k_OdfBinsPhi1/PHI/Phi2) documenting the fixed 5-degree
Bunge-Euler MATLAB ODF grid layout (186624 bins).

FIX 3: Guard in simulateMTR that pgrf_result.mtrIndex.size() == N; throws
std::runtime_error if the PGRF result dimensions are inconsistent.

FIX 4: CSV loop now uses sim.nx*sim.ny*sim.nz to tie the bound to the
actual result rather than the pre-call N.

FIX 5: Add phi2 range check and phi1/phi/phi2 size-equality checks to the
statistical driver test.

FIX 6: Remove unused <numbers> include from main.cpp; ODFSampler.hpp
retained (ODFComponent is used directly).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: imikejackson

src/LibMTRSim/MTRSimDriver.cpp

@@ -85,6 +85,10 @@ MTRSimResult simulateMTR(const SimulationParams& params,
   PGRFSimulation pgrf{rng};
   const PGRFResult pgrf_result = pgrf.run(params); // throws on bad dims
 
+  if (static_cast<int>(pgrf_result.mtrIndex.size()) != N) {
+    throw std::runtime_error("simulateMTR: PGRF result size does not match grid dimensions");
+  }
+
   // 2. Sample N orientations per component against the uniform reference.
   const ODFComponent uniformOdf = buildUniformODF(n1, nPHI, n2);
   const int numComponents = static_cast<int>(odfComponents.size());

src/app/main.cpp

@@ -14,7 +14,6 @@
 
 #include <cmath>
 #include <fstream>
-#include <numbers>
 #include <random>
 #include <stdexcept>
 #include <string>
@@ -203,19 +202,20 @@ int main(int argc, char **argv) {
   spdlog::info("  dx={}  dy={}  dz={} [mm]", params.dx, params.dy, params.dz);
 
   // ── Build spatial coordinate matrix
-  // ────────────────────────────────────────── Ordering matches
-  // simulate_MTRs.m: z outer → x middle → y inner.
-  //   s(k) = [j*dx,  i*dy,  zix*dz]
-  //          j ∈ [1,nx],  i ∈ [1,ny],  zix ∈ [1,nz]
+  // ──────────────────────────────────────────
+  // Row order matches simulateMTR's SIMPLNX z,y,x output: z slowest, x fastest.
+  //   k = iz*(ny*nx) + iy*nx + ix
+  //   spatialCoords(k) = [(ix+1)*dx, (iy+1)*dy, (iz+1)*dz]
+  // 1-based coordinate values are preserved (matching the original convention).
   Eigen::MatrixXd spatialCoords(N, 3);
   {
-    int k = 0;
-    for (int zix = 1; zix <= nz; ++zix) {
-      for (int j = 1; j <= nx; ++j) {
-        for (int i = 1; i <= ny; ++i, ++k) {
-          spatialCoords(k, 0) = j * params.dx;
-          spatialCoords(k, 1) = i * params.dy;
-          spatialCoords(k, 2) = zix * params.dz;
+    for (int iz = 0; iz < nz; ++iz) {
+      for (int iy = 0; iy < ny; ++iy) {
+        for (int ix = 0; ix < nx; ++ix) {
+          const int k = iz * (ny * nx) + iy * nx + ix;
+          spatialCoords(k, 0) = (ix + 1) * params.dx;
+          spatialCoords(k, 1) = (iy + 1) * params.dy;
+          spatialCoords(k, 2) = (iz + 1) * params.dz;
         }
       }
     }
@@ -239,8 +239,14 @@ int main(int argc, char **argv) {
 
   // ── Run full MTR simulation (PGRF + ODF sampling + remap to z,y,x order)
   // ─────────────────────────────────────────────────────────────────────────
+  // The MATLAB ODF HDF5 layout is a fixed 5-degree Bunge-Euler grid:
+  // 72 (phi1) x 36 (PHI) x 72 (phi2) = 186624 bins.
+  constexpr int k_OdfBinsPhi1 = 72;
+  constexpr int k_OdfBinsPHI = 36;
+  constexpr int k_OdfBinsPhi2 = 72;
+
   spdlog::info("Running MTR simulation...");
-  mtrsim::MTRSimResult sim = mtrsim::simulateMTR(params, odfComponents, rng, 72, 36, 72);
+  mtrsim::MTRSimResult sim = mtrsim::simulateMTR(params, odfComponents, rng, k_OdfBinsPhi1, k_OdfBinsPHI, k_OdfBinsPhi2);
   spdlog::info("MTR simulation complete.");
 
   Eigen::VectorXd phi1Vec = Eigen::Map<Eigen::VectorXd>(sim.phi1.data(), static_cast<Eigen::Index>(sim.phi1.size()));
@@ -275,7 +281,9 @@ int main(int argc, char **argv) {
     csv << std::fixed;
     csv.precision(6);
 
-    for (int i = 0; i < N; ++i) {
+    // Use sim dimensions to tie loop bounds to the actual result.
+    const int simN = sim.nx * sim.ny * sim.nz;
+    for (int i = 0; i < simN; ++i) {
       csv << spatialCoords(i, 0) << ',' << spatialCoords(i, 1) << ','
           << spatialCoords(i, 2) << ',' << phi1Vec[i] << ',' << phiVec[i] << ','
           << phi2Vec[i] << ',' << sim.mtrIndex[static_cast<std::size_t>(i)] << '\n';

tests/test_mtrsim_driver.cpp

@@ -88,6 +88,11 @@ TEST_CASE("simulateMTR reproduces target volume fractions (statistical)", "[mtrs
   REQUIRE(static_cast<double>(counts[1]) / N == Approx(0.35).margin(0.05));
   REQUIRE(static_cast<double>(counts[2]) / N == Approx(0.35).margin(0.05));
 
+  REQUIRE(static_cast<int>(r.phi1.size()) == N);
+  REQUIRE(static_cast<int>(r.phi.size()) == N);
+  REQUIRE(static_cast<int>(r.phi2.size()) == N);
+
   for (double a : r.phi1) { REQUIRE(a >= 0.0); REQUIRE(a <= 2.0 * std::numbers::pi); }
   for (double a : r.phi) { REQUIRE(a >= 0.0); REQUIRE(a <= std::numbers::pi); }
+  for (double a : r.phi2) { REQUIRE(a >= 0.0); REQUIRE(a <= 2.0 * std::numbers::pi); }
 }