@@ -124,86 +124,66 @@ inline void BuildOrientationTestData(DataStructure& ds, const ShapeType& cellSha
124124 auto * phasesArray = DataArray<int32>::Create (ds, " Phases" , phasesDataStore, amId);
125125 auto & phasesStore = phasesArray->getDataStoreRef ();
126126
127- // Quaternion Hamilton product: result = a * b, where q = (w, x, y, z)
128- auto quatMul = [](const std::array<float32, 4 >& a, const std::array<float32, 4 >& b) -> std::array<float32, 4 > {
129- return {a[0 ] * b[0 ] - a[1 ] * b[1 ] - a[2 ] * b[2 ] - a[3 ] * b[3 ], a[0 ] * b[1 ] + a[1 ] * b[0 ] + a[2 ] * b[3 ] - a[3 ] * b[2 ],
130- a[0 ] * b[2 ] - a[1 ] * b[3 ] + a[2 ] * b[0 ] + a[3 ] * b[1 ], a[0 ] * b[3 ] + a[1 ] * b[2 ] - a[2 ] * b[1 ] + a[3 ] * b[0 ]};
131- };
132-
133- constexpr float32 k_Pi = 3 .14159265358979323846f ;
127+ constexpr float32 k_DegToRad = 3 .14159265358979323846f / 180 .0f ;
134128
135129 const usize blocksPerX = (dimX + blockSize - 1 ) / blockSize;
136130 const usize blocksPerY = (dimY + blockSize - 1 ) / blockSize;
137131 const usize blocksPerZ = (dimZ + blockSize - 1 ) / blockSize;
138132 const usize numBlocks = blocksPerX * blocksPerY * blocksPerZ;
139133
140- // Pre-compute one quaternion per block. Two strategies are needed because
141- // EBSD compares full misorientations (24 cubic symmetry operators fold
142- // orientation space) while CAxis compares only C-axis directions (a 2D
143- // quantity on the hemisphere).
134+ // Quaternion Hamilton product: result = a * b, where q = (w, x, y, z)
135+ auto quatMul = [](const std::array<float32, 4 >& a, const std::array<float32, 4 >& b) -> std::array<float32, 4 > {
136+ return {a[0 ] * b[0 ] - a[1 ] * b[1 ] - a[2 ] * b[2 ] - a[3 ] * b[3 ], a[0 ] * b[1 ] + a[1 ] * b[0 ] + a[2 ] * b[3 ] - a[3 ] * b[2 ],
137+ a[0 ] * b[2 ] - a[1 ] * b[3 ] + a[2 ] * b[0 ] + a[3 ] * b[1 ], a[0 ] * b[3 ] + a[1 ] * b[2 ] - a[2 ] * b[1 ] + a[3 ] * b[0 ]};
138+ };
139+
144140 std::vector<std::array<float32, 4 >> blockQuats (numBlocks);
145141
146- if (crystalStructure == 0 ) // Hexagonal_High → CAxis comparison uses C-axis directions
142+ // Z-layer orientation scheme (shared by EBSD and CAxis):
143+ // All blocks in the same Z-layer share a single X-axis rotation angle.
144+ // This produces 3 horizontal layers of identical orientations:
145+ // z=0: 0° rotation → q = [1, 0, 0, 0] c-axis = [0, 0, 1]
146+ // z=1: 30° rotation → q = [0.966, 0.259, 0, 0] c-axis = [0, 0.5, 0.866]
147+ // z=2: 60° rotation → q = [0.866, 0.5, 0, 0] c-axis = [0, 0.866, 0.5]
148+ //
149+ // Adjacent layers differ by 30°, well above the 5° tolerance → no merge.
150+ // Within each layer, all blocks share the same angle → they merge.
151+ //
152+ // Merge pair override (non-periodic only):
153+ // Block (1,1,1) at center of z=1 is set to 0° instead of 30°.
154+ // It merges with its z=0 neighbor (1,1,0) while staying separate
155+ // from the other z=1 blocks (30° difference → no merge).
156+ //
157+ // Expected features (3x3x3 grid):
158+ // Base: 3 (z=0 layer + center pillar, z=1 minus pillar, z=2 layer)
159+ // Periodic: 2 (z=0 and z=2 share 0° via wrapping → merge, z=1 separate)
160+ constexpr float32 k_LayerAngles[] = {0 .0f , 30 .0f , 60 .0f };
161+
162+ for (usize bz = 0 ; bz < blocksPerZ; bz++)
147163 {
148- // Fibonacci hemisphere: distributes block C-axes uniformly on the upper
149- // hemisphere with ~8-9 degree minimum separation for ≤125 blocks, which
150- // comfortably exceeds the 5-degree tolerance used in tests.
151- //
152- // The CAxis filter computes the sample-frame C-axis as:
153- // c_sample = oMatrix.transpose() * [0,0,1]
154- // where oMatrix = rotationMatrix(q). This equals q^{-1} * [0,0,1].
155- // So we need q^{-1} * [0,0,1] = fibonacci_point, meaning
156- // q = Ry(-theta) * Rz(-phi)
157- // (the inverse of the rotation FROM [0,0,1] TO the hemisphere point).
158- constexpr float32 k_GoldenAngle = 2 .399963229728653f ; // pi * (3 - sqrt(5))
159- for (usize i = 0 ; i < numBlocks; i++)
160- {
161- const float32 cosTheta = 1 .0f - (static_cast <float32>(i) + 0 .5f ) / static_cast <float32>(numBlocks);
162- const float32 theta = std::acos (std::clamp (cosTheta, 0 .0f , 1 .0f ));
163- const float32 phi = static_cast <float32>(i) * k_GoldenAngle;
164-
165- const float32 halfPhi = phi * 0 .5f ;
166- const float32 halfTheta = theta * 0 .5f ;
167- const std::array<float32, 4 > qyNeg = {std::cos (halfTheta), 0 .0f , -std::sin (halfTheta), 0 .0f };
168- const std::array<float32, 4 > qzNeg = {std::cos (halfPhi), 0 .0f , 0 .0f , -std::sin (halfPhi)};
169- blockQuats[i] = quatMul (qyNeg, qzNeg);
170- }
171- }
172- else // Cubic_High → EBSD comparison uses full misorientation
173- {
174- // Composed rotations around X, Y, and (1,1,0)/sqrt(2). Each block
175- // index contributes an independent 14-degree step. Adjacent blocks
176- // differ by 14 degrees in one component — well above the 5-degree
177- // tolerance under cubic symmetry's 24-operator reduction.
178- // NOTE: Step must be chosen so that (blocksPerAxis-1)*step != 90,
179- // because 90-degree rotations around <100> are cubic symmetry
180- // operators. With periodic wrapping (7 effective blocks), 6*14=84
181- // has misorientation |90-84|=6 degrees > 5-degree tolerance.
182- constexpr float32 k_Step = 14 .0f * (k_Pi / 180 .0f );
183- constexpr float32 k_InvSqrt2 = 0 .70710678118654752f ;
184-
185- for (usize bz = 0 ; bz < blocksPerZ; bz++)
164+ const usize layerIdx = std::min (bz, static_cast <usize>(2 ));
165+ const float32 halfAngle = k_LayerAngles[layerIdx] * k_DegToRad * 0 .5f ;
166+ // EBSDlib quaternion layout: (x, y, z, w) — Vector-Scalar order
167+ const std::array<float32, 4 > layerQuat = {std::sin (halfAngle), 0 .0f , 0 .0f , std::cos (halfAngle)};
168+
169+ for (usize by = 0 ; by < blocksPerY; by++)
186170 {
187- for (usize by = 0 ; by < blocksPerY; by ++)
171+ for (usize bx = 0 ; bx < blocksPerX; bx ++)
188172 {
189- for (usize bx = 0 ; bx < blocksPerX; bx++)
190- {
191- const float32 ax = static_cast <float32>(bx) * k_Step;
192- const float32 ay = static_cast <float32>(by) * k_Step;
193- const float32 az = static_cast <float32>(bz) * k_Step;
194-
195- const std::array<float32, 4 > qx = {std::cos (ax * 0 .5f ), std::sin (ax * 0 .5f ), 0 .0f , 0 .0f };
196- const std::array<float32, 4 > qyRot = {std::cos (ay * 0 .5f ), 0 .0f , std::sin (ay * 0 .5f ), 0 .0f };
197- const float32 halfAz = az * 0 .5f ;
198- const std::array<float32, 4 > qd = {std::cos (halfAz), std::sin (halfAz) * k_InvSqrt2, std::sin (halfAz) * k_InvSqrt2, 0 .0f };
199-
200- const usize blockIdx = bz * blocksPerY * blocksPerX + by * blocksPerX + bx;
201- blockQuats[blockIdx] = quatMul (qd, quatMul (qyRot, qx));
202- }
173+ const usize blockIdx = bz * blocksPerY * blocksPerX + by * blocksPerX + bx;
174+ blockQuats[blockIdx] = layerQuat;
203175 }
204176 }
205177 }
206178
179+ // Merge pair: block (1,1,1) gets z=0 angle (0°) instead of z=1 angle (30°).
180+ // It merges downward into the z=0 layer through face neighbor (1,1,0).
181+ if (!wrapBoundary && blocksPerX >= 3 && blocksPerY >= 3 && blocksPerZ >= 3 )
182+ {
183+ const usize idx_111 = 1 * blocksPerY * blocksPerX + 1 * blocksPerX + 1 ;
184+ blockQuats[idx_111] = blockQuats[0 ]; // Set to 0° (z=0 layer angle)
185+ }
186+
207187 for (usize z = 0 ; z < dimZ; z++)
208188 {
209189 for (usize y = 0 ; y < dimY; y++)
0 commit comments