|
3 | 3 | #include "simplnx/Common/Numbers.hpp" |
4 | 4 | #include "simplnx/DataStructure/DataGroup.hpp" |
5 | 5 | #include "simplnx/DataStructure/Geometry/IGridGeometry.hpp" |
| 6 | +#include "simplnx/Utilities/AlgorithmDispatch.hpp" |
6 | 7 | #include "simplnx/Utilities/FilterUtilities.hpp" |
7 | 8 | #include "simplnx/Utilities/MaskCompareUtilities.hpp" |
8 | 9 |
|
@@ -41,6 +42,15 @@ Result<> AlignSectionsMisorientation::operator()() |
41 | 42 | // ----------------------------------------------------------------------------- |
42 | 43 | Result<> AlignSectionsMisorientation::findShifts(std::vector<int64_t>& xShifts, std::vector<int64_t>& yShifts) |
43 | 44 | { |
| 45 | + { |
| 46 | + const auto& quatsCheck = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->quatsArrayPath); |
| 47 | + const auto& cellPhasesCheck = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->cellPhasesArrayPath); |
| 48 | + if(ForceOocAlgorithm() || IsOutOfCore(quatsCheck) || IsOutOfCore(cellPhasesCheck)) |
| 49 | + { |
| 50 | + return findShiftsOoc(xShifts, yShifts); |
| 51 | + } |
| 52 | + } |
| 53 | + |
44 | 54 | std::unique_ptr<MaskCompareUtilities::MaskCompare> maskCompare = nullptr; |
45 | 55 | if(m_InputValues->UseMask) |
46 | 56 | { |
@@ -298,3 +308,217 @@ Result<> AlignSectionsMisorientation::findShifts(std::vector<int64_t>& xShifts, |
298 | 308 |
|
299 | 309 | return {}; |
300 | 310 | } |
| 311 | + |
| 312 | +// ----------------------------------------------------------------------------- |
| 313 | +// OOC-optimized findShifts: buffers 2 adjacent Z-slices of quats, cellPhases, |
| 314 | +// and mask into local vectors before the convergence loop, eliminating random |
| 315 | +// chunk-based DataStore access. |
| 316 | +// ----------------------------------------------------------------------------- |
| 317 | +Result<> AlignSectionsMisorientation::findShiftsOoc(std::vector<int64_t>& xShifts, std::vector<int64_t>& yShifts) |
| 318 | +{ |
| 319 | + std::unique_ptr<MaskCompareUtilities::MaskCompare> maskCompare = nullptr; |
| 320 | + if(m_InputValues->UseMask) |
| 321 | + { |
| 322 | + try |
| 323 | + { |
| 324 | + maskCompare = MaskCompareUtilities::InstantiateMaskCompare(m_DataStructure, m_InputValues->MaskArrayPath); |
| 325 | + } catch(const std::out_of_range& exception) |
| 326 | + { |
| 327 | + std::string message = fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", m_InputValues->MaskArrayPath.toString()); |
| 328 | + return MakeErrorResult(-53900, message); |
| 329 | + } |
| 330 | + } |
| 331 | + |
| 332 | + auto* gridGeom = m_DataStructure.getDataAs<IGridGeometry>(m_InputValues->ImageGeometryPath); |
| 333 | + |
| 334 | + const auto& cellPhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->cellPhasesArrayPath); |
| 335 | + const auto& quats = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->quatsArrayPath); |
| 336 | + const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->crystalStructuresArrayPath); |
| 337 | + auto& cellPhasesStore = cellPhases.getDataStoreRef(); |
| 338 | + auto& quatsStore = quats.getDataStoreRef(); |
| 339 | + |
| 340 | + SizeVec3 udims = gridGeom->getDimensions(); |
| 341 | + |
| 342 | + std::array<int64_t, 3> dims = { |
| 343 | + static_cast<int64_t>(udims[0]), |
| 344 | + static_cast<int64_t>(udims[1]), |
| 345 | + static_cast<int64_t>(udims[2]), |
| 346 | + }; |
| 347 | + |
| 348 | + std::vector<ebsdlib::LaueOps::Pointer> orientationOps = ebsdlib::LaueOps::GetAllOrientationOps(); |
| 349 | + |
| 350 | + std::vector<bool> misorients(dims[0] * dims[1], false); |
| 351 | + |
| 352 | + const auto halfDim0 = static_cast<int64_t>(dims[0] * 0.5f); |
| 353 | + const auto halfDim1 = static_cast<int64_t>(dims[1] * 0.5f); |
| 354 | + |
| 355 | + double deg2Rad = (nx::core::numbers::pi / 180.0); |
| 356 | + ThrottledMessenger throttledMessenger = getMessageHelper().createThrottledMessenger(); |
| 357 | + |
| 358 | + const int64_t sliceVoxels = dims[0] * dims[1]; |
| 359 | + |
| 360 | + // Buffers for 2 Z-slices: reference (slice+1) and current (slice) |
| 361 | + std::vector<float32> refQuatsBuf(sliceVoxels * 4); |
| 362 | + std::vector<float32> curQuatsBuf(sliceVoxels * 4); |
| 363 | + std::vector<int32_t> refPhasesBuf(sliceVoxels); |
| 364 | + std::vector<int32_t> curPhasesBuf(sliceVoxels); |
| 365 | + std::vector<uint8_t> refMaskBuf; |
| 366 | + std::vector<uint8_t> curMaskBuf; |
| 367 | + if(m_InputValues->UseMask) |
| 368 | + { |
| 369 | + refMaskBuf.resize(sliceVoxels, 1); |
| 370 | + curMaskBuf.resize(sliceVoxels, 1); |
| 371 | + } |
| 372 | + |
| 373 | + // Optional output stores |
| 374 | + AbstractDataStore<uint32>* slicesStorePtr = nullptr; |
| 375 | + AbstractDataStore<int64>* relativeShiftsStorePtr = nullptr; |
| 376 | + AbstractDataStore<int64>* cumulativeShiftsStorePtr = nullptr; |
| 377 | + if(m_InputValues->StoreAlignmentShifts) |
| 378 | + { |
| 379 | + slicesStorePtr = &m_DataStructure.getDataAs<UInt32Array>(m_InputValues->SlicesArrayPath)->getDataStoreRef(); |
| 380 | + relativeShiftsStorePtr = &m_DataStructure.getDataAs<Int64Array>(m_InputValues->RelativeShiftsArrayPath)->getDataStoreRef(); |
| 381 | + cumulativeShiftsStorePtr = &m_DataStructure.getDataAs<Int64Array>(m_InputValues->CumulativeShiftsArrayPath)->getDataStoreRef(); |
| 382 | + } |
| 383 | + |
| 384 | + for(int64_t iter = 1; iter < dims[2]; iter++) |
| 385 | + { |
| 386 | + if(m_ShouldCancel) |
| 387 | + { |
| 388 | + return {}; |
| 389 | + } |
| 390 | + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Determining Shifts || {:.2f}% Complete", CalculatePercentComplete(iter, dims[2])); }); |
| 391 | + if(getCancel()) |
| 392 | + { |
| 393 | + return {}; |
| 394 | + } |
| 395 | + |
| 396 | + int64_t slice = (dims[2] - 1) - iter; |
| 397 | + |
| 398 | + // Buffer reference slice (slice+1) and current slice (slice) |
| 399 | + int64_t refOffset = (slice + 1) * sliceVoxels; |
| 400 | + int64_t curOffset = slice * sliceVoxels; |
| 401 | + |
| 402 | + for(int64_t idx = 0; idx < sliceVoxels; idx++) |
| 403 | + { |
| 404 | + refPhasesBuf[idx] = cellPhasesStore[refOffset + idx]; |
| 405 | + curPhasesBuf[idx] = cellPhasesStore[curOffset + idx]; |
| 406 | + } |
| 407 | + for(int64_t idx = 0; idx < sliceVoxels * 4; idx++) |
| 408 | + { |
| 409 | + refQuatsBuf[idx] = quatsStore[refOffset * 4 + idx]; |
| 410 | + curQuatsBuf[idx] = quatsStore[curOffset * 4 + idx]; |
| 411 | + } |
| 412 | + if(m_InputValues->UseMask) |
| 413 | + { |
| 414 | + for(int64_t idx = 0; idx < sliceVoxels; idx++) |
| 415 | + { |
| 416 | + refMaskBuf[idx] = maskCompare->isTrue(refOffset + idx) ? 1 : 0; |
| 417 | + curMaskBuf[idx] = maskCompare->isTrue(curOffset + idx) ? 1 : 0; |
| 418 | + } |
| 419 | + } |
| 420 | + |
| 421 | + float minDisorientation = std::numeric_limits<float>::max(); |
| 422 | + int64_t oldxshift = -1; |
| 423 | + int64_t oldyshift = -1; |
| 424 | + int64_t newxshift = 0; |
| 425 | + int64_t newyshift = 0; |
| 426 | + |
| 427 | + std::fill(misorients.begin(), misorients.end(), false); |
| 428 | + |
| 429 | + float misorientationTolerance = static_cast<float>(m_InputValues->misorientationTolerance * deg2Rad); |
| 430 | + |
| 431 | + while(newxshift != oldxshift || newyshift != oldyshift) |
| 432 | + { |
| 433 | + oldxshift = newxshift; |
| 434 | + oldyshift = newyshift; |
| 435 | + for(int32_t j = -3; j < 4; j++) |
| 436 | + { |
| 437 | + for(int32_t k = -3; k < 4; k++) |
| 438 | + { |
| 439 | + float disorientation = 0.0f; |
| 440 | + float count = 0.0f; |
| 441 | + int64_t xIdx = k + oldxshift + halfDim0; |
| 442 | + int64_t yIdx = j + oldyshift + halfDim1; |
| 443 | + int64_t idx = (dims[0] * yIdx) + xIdx; |
| 444 | + if(!misorients[idx] && llabs(k + oldxshift) < halfDim0 && llabs(j + oldyshift) < halfDim1) |
| 445 | + { |
| 446 | + for(int64_t l = 0; l < dims[1]; l = l + 4) |
| 447 | + { |
| 448 | + for(int64_t n = 0; n < dims[0]; n = n + 4) |
| 449 | + { |
| 450 | + if((l + j + oldyshift) >= 0 && (l + j + oldyshift) < dims[1] && (n + k + oldxshift) >= 0 && (n + k + oldxshift) < dims[0]) |
| 451 | + { |
| 452 | + count++; |
| 453 | + // Local buffer indices (within-slice) |
| 454 | + int64_t refLocalIdx = l * dims[0] + n; |
| 455 | + int64_t curLocalIdx = (l + j + oldyshift) * dims[0] + (n + k + oldxshift); |
| 456 | + |
| 457 | + bool maskOk = !m_InputValues->UseMask || (refMaskBuf[refLocalIdx] != 0 && curMaskBuf[curLocalIdx] != 0); |
| 458 | + if(maskOk) |
| 459 | + { |
| 460 | + float angle = std::numeric_limits<float>::max(); |
| 461 | + if(refPhasesBuf[refLocalIdx] > 0 && curPhasesBuf[curLocalIdx] > 0) |
| 462 | + { |
| 463 | + ebsdlib::QuatD quat1(refQuatsBuf[refLocalIdx * 4], refQuatsBuf[refLocalIdx * 4 + 1], refQuatsBuf[refLocalIdx * 4 + 2], refQuatsBuf[refLocalIdx * 4 + 3]); |
| 464 | + auto laueClass1 = static_cast<int32_t>(crystalStructures[refPhasesBuf[refLocalIdx]]); |
| 465 | + ebsdlib::QuatD quat2(curQuatsBuf[curLocalIdx * 4], curQuatsBuf[curLocalIdx * 4 + 1], curQuatsBuf[curLocalIdx * 4 + 2], curQuatsBuf[curLocalIdx * 4 + 3]); |
| 466 | + auto laueClass2 = static_cast<int32_t>(crystalStructures[curPhasesBuf[curLocalIdx]]); |
| 467 | + if(laueClass1 == laueClass2 && laueClass1 < static_cast<uint32_t>(orientationOps.size())) |
| 468 | + { |
| 469 | + ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2); |
| 470 | + angle = axisAngle[3]; |
| 471 | + } |
| 472 | + } |
| 473 | + if(angle > misorientationTolerance) |
| 474 | + { |
| 475 | + disorientation++; |
| 476 | + } |
| 477 | + } |
| 478 | + if(m_InputValues->UseMask) |
| 479 | + { |
| 480 | + if(refMaskBuf[refLocalIdx] != 0 && curMaskBuf[curLocalIdx] == 0) |
| 481 | + { |
| 482 | + disorientation++; |
| 483 | + } |
| 484 | + if(refMaskBuf[refLocalIdx] == 0 && curMaskBuf[curLocalIdx] != 0) |
| 485 | + { |
| 486 | + disorientation++; |
| 487 | + } |
| 488 | + } |
| 489 | + } |
| 490 | + } |
| 491 | + } |
| 492 | + disorientation = disorientation / count; |
| 493 | + xIdx = k + oldxshift + halfDim0; |
| 494 | + yIdx = j + oldyshift + halfDim1; |
| 495 | + idx = (dims[0] * yIdx) + xIdx; |
| 496 | + misorients[idx] = true; |
| 497 | + if(disorientation < minDisorientation || (disorientation == minDisorientation && ((llabs(k + oldxshift) < llabs(newxshift)) || (llabs(j + oldyshift) < llabs(newyshift))))) |
| 498 | + { |
| 499 | + newxshift = k + oldxshift; |
| 500 | + newyshift = j + oldyshift; |
| 501 | + minDisorientation = disorientation; |
| 502 | + } |
| 503 | + } |
| 504 | + } |
| 505 | + } |
| 506 | + } |
| 507 | + xShifts[iter] = xShifts[iter - 1] + newxshift; |
| 508 | + yShifts[iter] = yShifts[iter - 1] + newyshift; |
| 509 | + |
| 510 | + if(m_InputValues->StoreAlignmentShifts) |
| 511 | + { |
| 512 | + usize xIndex = iter * 2; |
| 513 | + usize yIndex = (iter * 2) + 1; |
| 514 | + (*slicesStorePtr)[xIndex] = slice; |
| 515 | + (*slicesStorePtr)[yIndex] = slice + 1; |
| 516 | + (*relativeShiftsStorePtr)[xIndex] = newxshift; |
| 517 | + (*relativeShiftsStorePtr)[yIndex] = newyshift; |
| 518 | + (*cumulativeShiftsStorePtr)[xIndex] = xShifts[iter]; |
| 519 | + (*cumulativeShiftsStorePtr)[yIndex] = yShifts[iter]; |
| 520 | + } |
| 521 | + } |
| 522 | + |
| 523 | + return {}; |
| 524 | +} |
0 commit comments