Skip to content

Commit b107269

Browse files
imikejacksonclaude
andcommitted
feat(filter): implement MTRSim algorithm execute + statistical wiring test
Wire MTRSimFilter to the LibMTRSim simulateMTR entry point: reconstruct ODFComponents from the selected Float64 cell arrays, build SimulationParams from filter inputs, run the simulation, and write MTRIds + Eulers into the output geometry cell AttributeMatrix. Polar coloring (Task 8) is deferred. Add preflight non-negativity/range hardening (error -13007) requiring each Volume Fraction value to lie in [0, 1]. Add an end-to-end execute test (100x100 = 10000 voxels, seed 42, 3 components) that verifies array contract (types/components/tuples), id set {1,2,3}, finite Bunge-bounded Eulers, recorded seed, and loose volume fractions. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
1 parent 38bb3e0 commit b107269

3 files changed

Lines changed: 198 additions & 3 deletions

File tree

src/MTRSim/Filters/Algorithms/MTRSim.cpp

Lines changed: 88 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,17 @@
11
#include "MTRSim.hpp"
22

3+
#include "simplnx/DataStructure/DataArray.hpp"
4+
#include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
5+
6+
#include "LibMTRSim/MTRSimDriver.hpp"
7+
8+
#include <fmt/format.h>
9+
10+
#include <cmath>
11+
#include <exception>
12+
#include <random>
13+
#include <vector>
14+
315
using namespace nx::core;
416

517
// -----------------------------------------------------------------------------
@@ -17,8 +29,81 @@ MTRSim::~MTRSim() noexcept = default;
1729
// -----------------------------------------------------------------------------
1830
Result<> MTRSim::operator()()
1931
{
20-
// NOTE: The algorithm body is implemented in a later task. For now this is a
21-
// no-op stub so the filter scaffolding (preflight + parameter validation) can
22-
// be built and tested independently.
32+
m_MessageHandler(IFilter::Message::Type::Info, "Reading ODF geometry...");
33+
34+
// 1. ODF grid geometry (degrees) from the input ODF ImageGeom.
35+
const auto& odfGeom = m_DataStructure.getDataRefAs<ImageGeom>(m_InputValues->inputOdfGeometryPath);
36+
const SizeVec3 odfDims = odfGeom.getDimensions(); // X=phi2, Y=PHI, Z=phi1
37+
const FloatVec3 odfSpacing = odfGeom.getSpacing(); // degrees
38+
const int n2 = static_cast<int>(odfDims[0]);
39+
const int nPHI = static_cast<int>(odfDims[1]);
40+
const int n1 = static_cast<int>(odfDims[2]);
41+
42+
// 2. Reconstruct ODFComponents from the selected Float64 cell arrays. The
43+
// store is row-major ZYX (phi1 slowest, phi2 fastest), exactly the flat
44+
// layout gridToODFComponent expects, so we copy values in their natural order.
45+
std::vector<mtrsim::ODFComponent> components;
46+
components.reserve(m_InputValues->odfComponentPaths.size());
47+
for(const auto& path : m_InputValues->odfComponentPaths)
48+
{
49+
const auto& arr = m_DataStructure.getDataRefAs<Float64Array>(path);
50+
const auto& store = arr.getDataStoreRef();
51+
std::vector<double> values(store.begin(), store.end());
52+
components.push_back(mtrsim::gridToODFComponent(values, n1, nPHI, n2, static_cast<double>(odfSpacing[2]), static_cast<double>(odfSpacing[1]), static_cast<double>(odfSpacing[0])));
53+
}
54+
55+
// 3. SimulationParams from filter inputs.
56+
mtrsim::SimulationParams params;
57+
params.xLen = m_InputValues->physicalSize[0];
58+
params.yLen = m_InputValues->physicalSize[1];
59+
params.zLen = m_InputValues->physicalSize[2];
60+
params.dx = m_InputValues->physicalSpacing[0];
61+
params.dy = m_InputValues->physicalSpacing[1];
62+
params.dz = m_InputValues->physicalSpacing[2];
63+
params.volumeFractions = m_InputValues->volumeFractions[0]; // 1 row
64+
params.thetaList = m_InputValues->thetaList;
65+
params.seed = m_InputValues->seed;
66+
67+
if(m_ShouldCancel)
68+
{
69+
return {};
70+
}
71+
72+
// 4. Run simulation (SIMPLNX z,y,x order out).
73+
m_MessageHandler(IFilter::Message::Type::Info, "Running MTR simulation (this may take a while for large volumes)...");
74+
std::mt19937_64 rng(m_InputValues->seed);
75+
mtrsim::MTRSimResult sim;
76+
try
77+
{
78+
sim = mtrsim::simulateMTR(params, components, rng, n1, nPHI, n2);
79+
} catch(const std::exception& e)
80+
{
81+
return MakeErrorResult(-13050, fmt::format("MTR simulation failed: {}", e.what()));
82+
}
83+
84+
if(m_ShouldCancel)
85+
{
86+
return {};
87+
}
88+
89+
// 5. Write MTR ids + Euler arrays (output geometry cell AM).
90+
const DataPath cellAm = m_InputValues->outputGeometryPath.createChildPath(m_InputValues->cellAttrMatName);
91+
auto& mtrIds = m_DataStructure.getDataRefAs<Int32Array>(cellAm.createChildPath(m_InputValues->mtrIdsArrayName));
92+
auto& eulers = m_DataStructure.getDataRefAs<Float32Array>(cellAm.createChildPath(m_InputValues->eulersArrayName));
93+
auto& mtrStore = mtrIds.getDataStoreRef();
94+
auto& eulerStore = eulers.getDataStoreRef();
95+
96+
const std::size_t N = sim.mtrIndex.size();
97+
for(std::size_t i = 0; i < N; ++i)
98+
{
99+
mtrStore[i] = sim.mtrIndex[i];
100+
eulerStore[i * 3 + 0] = static_cast<float>(sim.phi1[i]);
101+
eulerStore[i * 3 + 1] = static_cast<float>(sim.phi[i]);
102+
eulerStore[i * 3 + 2] = static_cast<float>(sim.phi2[i]);
103+
}
104+
105+
m_MessageHandler(IFilter::Message::Type::Info, "MTR simulation complete.");
106+
// Polar coloring (Task 8) is handled in a later task; the optional array (if
107+
// requested) is left created-but-unfilled here and populated next task.
23108
return {};
24109
}

src/MTRSim/Filters/MTRSimFilter.cpp

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -161,6 +161,13 @@ IFilter::PreflightResult MTRSimFilter::preflightImpl(const DataStructure& dataSt
161161
{
162162
return {MakeErrorResult<OutputActions>(-13003, fmt::format("Volume Fraction values must sum to 1.0 (got {:.4f}).", vfSum))};
163163
}
164+
for(double v : pVolumeFractions[0])
165+
{
166+
if(v < 0.0 || v > 1.0)
167+
{
168+
return {MakeErrorResult<OutputActions>(-13007, fmt::format("Each Volume Fraction value must be in the range [0, 1] (got {:.4f}).", v))};
169+
}
170+
}
164171
if(pThetaList.size() < numComponents - 1)
165172
{
166173
return {MakeErrorResult<OutputActions>(-13004, fmt::format("Theta List needs at least {} rows (components - 1).", numComponents - 1))};

test/MTRSimTest.cpp

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@
2323

2424
#include <fmt/format.h>
2525

26+
#include <array>
27+
#include <cmath>
2628
#include <vector>
2729

2830
using namespace nx::core;
@@ -116,6 +118,107 @@ TEST_CASE("MTRSim::MTRSimFilter: Rejects mismatched Volume Fraction column count
116118
SIMPLNX_RESULT_REQUIRE_INVALID(preflightResult.outputActions);
117119
}
118120

121+
TEST_CASE("MTRSim::MTRSimFilter: Execute wires simulation to output arrays", "[MTRSim][MTRSimFilter]")
122+
{
123+
UnitTest::LoadPlugins();
124+
125+
DataStructure dataStructure;
126+
const std::vector<DataPath> compPaths = BuildOdfDataStructure(dataStructure, 3);
127+
// Fill every component array with a uniform value so ODF sampling is well-defined.
128+
for(const auto& path : compPaths)
129+
{
130+
auto& arr = dataStructure.getDataRefAs<Float64Array>(path);
131+
arr.fill(1.0);
132+
}
133+
134+
MTRSimFilter filter;
135+
Arguments args = MakeValidArgs(compPaths);
136+
// Modest domain: 100x100 = 10000 voxels (Z=0 -> single layer).
137+
args.insertOrAssign(MTRSimFilter::k_PhysicalSize_Key, std::vector<float32>{2.0f, 2.0f, 0.0f});
138+
args.insertOrAssign(MTRSimFilter::k_PhysicalSpacing_Key, std::vector<float32>{0.02f, 0.02f, 0.02f});
139+
args.insertOrAssign(MTRSimFilter::k_VolumeFractions_Key, DynamicTableParameter::ValueType{{0.30, 0.35, 0.35}});
140+
args.insertOrAssign(MTRSimFilter::k_UseSeed_Key, true);
141+
args.insertOrAssign(MTRSimFilter::k_SeedValue_Key, static_cast<uint64>(42));
142+
143+
auto preflightResult = filter.preflight(dataStructure, args);
144+
SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions);
145+
146+
auto executeResult = filter.execute(dataStructure, args);
147+
SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result);
148+
149+
// Output ImageGeom exists.
150+
const DataPath outGeomPath({"MTR Microstructure"});
151+
REQUIRE_NOTHROW(dataStructure.getDataRefAs<ImageGeom>(outGeomPath));
152+
153+
const DataPath cellAm = outGeomPath.createChildPath("Cell Data");
154+
const usize expectedTuples = 100 * 100;
155+
156+
// MTRIds: Int32, 1 component, 10000 tuples.
157+
auto& mtrIds = dataStructure.getDataRefAs<Int32Array>(cellAm.createChildPath("MTRIds"));
158+
REQUIRE(mtrIds.getNumberOfComponents() == 1);
159+
REQUIRE(mtrIds.getNumberOfTuples() == expectedTuples);
160+
161+
// Eulers: Float32, 3 components, 10000 tuples.
162+
auto& eulers = dataStructure.getDataRefAs<Float32Array>(cellAm.createChildPath("Eulers"));
163+
REQUIRE(eulers.getNumberOfComponents() == 3);
164+
REQUIRE(eulers.getNumberOfTuples() == expectedTuples);
165+
166+
// MTR ids in {1,2,3}; at least 2 distinct ids appear. Also accumulate empirical
167+
// volume fractions for a loose wiring check.
168+
const auto& mtrStore = mtrIds.getDataStoreRef();
169+
std::array<usize, 4> counts = {0, 0, 0, 0};
170+
for(usize i = 0; i < mtrStore.getSize(); ++i)
171+
{
172+
const int32 id = mtrStore[i];
173+
REQUIRE(id >= 1);
174+
REQUIRE(id <= 3);
175+
counts[static_cast<usize>(id)]++;
176+
}
177+
usize distinct = 0;
178+
for(usize id = 1; id <= 3; ++id)
179+
{
180+
if(counts[id] > 0)
181+
{
182+
distinct++;
183+
}
184+
}
185+
REQUIRE(distinct >= 2);
186+
187+
// Euler values finite and within Bunge bounds (interleaved 3/voxel).
188+
constexpr float twoPi = 2.0f * static_cast<float>(M_PI);
189+
constexpr float pi = static_cast<float>(M_PI);
190+
const auto& eulerStore = eulers.getDataStoreRef();
191+
for(usize t = 0; t < expectedTuples; ++t)
192+
{
193+
const float phi1 = eulerStore[t * 3 + 0];
194+
const float Phi = eulerStore[t * 3 + 1];
195+
const float phi2 = eulerStore[t * 3 + 2];
196+
REQUIRE(std::isfinite(phi1));
197+
REQUIRE(std::isfinite(Phi));
198+
REQUIRE(std::isfinite(phi2));
199+
REQUIRE(phi1 >= 0.0f);
200+
REQUIRE(phi1 <= twoPi);
201+
REQUIRE(Phi >= 0.0f);
202+
REQUIRE(Phi <= pi);
203+
REQUIRE(phi2 >= 0.0f);
204+
REQUIRE(phi2 <= twoPi);
205+
}
206+
207+
// Seed array records 42.
208+
auto& seedArray = dataStructure.getDataRefAs<UInt64Array>(DataPath({"MTRSim SeedValue"}));
209+
REQUIRE(seedArray[0] == 42);
210+
211+
// Loose volume-fraction wiring check (NOT a statistics check; rigorous VF
212+
// validation lives in the LibMTRSim statistical test). A 100x100 correlated
213+
// field has real variance, so use a generous margin of 0.12.
214+
const std::array<double, 4> targets = {0.0, 0.30, 0.35, 0.35};
215+
for(usize id = 1; id <= 3; ++id)
216+
{
217+
const double empirical = static_cast<double>(counts[id]) / static_cast<double>(expectedTuples);
218+
REQUIRE(empirical == Approx(targets[id]).margin(0.12));
219+
}
220+
}
221+
119222
TEST_CASE("MTRSim::MTRSimFilter: Rejects too few Theta List rows", "[MTRSim][MTRSimFilter][ErrorPath]")
120223
{
121224
UnitTest::LoadPlugins();

0 commit comments

Comments
 (0)