@@ -140,6 +140,8 @@ Image::ImageType::Pointer Image::read(const std::string& pathname) {
140140 img = orienter->GetOutput ();
141141 }
142142
143+ img = ImageUtils::make_axis_aligned (img);
144+
143145 return img;
144146}
145147
@@ -431,90 +433,6 @@ Image& Image::resize(Dims dims, Image::InterpolationType interp) {
431433 return resample (IdentityTransform::New (), origin (), dims, spacing, coordsys (), interp);
432434}
433435
434-
435- Image& Image::toAxisAligned (InterpolationType interp) {
436- // Check if image is already axis-aligned
437- auto direction = itk_image_->GetDirection ();
438- bool is_oblique = false ;
439- for (unsigned int i = 0 ; i < 3 ; i++) {
440- for (unsigned int j = 0 ; j < 3 ; j++) {
441- double expected = (i == j) ? 1.0 : 0.0 ;
442- if (std::abs (direction (i, j) - expected) > 1e-6 ) {
443- is_oblique = true ;
444- break ;
445- }
446- }
447- }
448-
449- // If already axis-aligned, return immediately
450- if (!is_oblique) {
451- return *this ;
452- }
453-
454- // The key insight: we need to create a transform that represents
455- // the DIFFERENCE between the oblique and axis-aligned coordinate systems
456-
457- using TransformType = itk::AffineTransform<double , 3 >;
458- auto transform = TransformType::New ();
459-
460- // The direction matrix rotates from index space to physical space
461- // To go from new axis-aligned indices to old oblique physical space,
462- // we need: newPhysical = identity * newIndex
463- // oldIndex = direction^-1 * oldPhysical
464- // So: oldIndex = direction^-1 * identity * newIndex = direction^-1 * newIndex
465-
466- TransformType::MatrixType matrix;
467- for (int i = 0 ; i < 3 ; i++) {
468- for (int j = 0 ; j < 3 ; j++) {
469- matrix (i, j) = direction (i, j);
470- }
471- }
472- transform->SetMatrix (matrix);
473-
474- // Get bounding box
475- auto size = itk_image_->GetLargestPossibleRegion ().GetSize ();
476- auto spacing = itk_image_->GetSpacing ();
477- auto origin = itk_image_->GetOrigin ();
478-
479- Point3 minPt, maxPt;
480- for (int i = 0 ; i < 3 ; i++) {
481- minPt[i] = std::numeric_limits<double >::max ();
482- maxPt[i] = std::numeric_limits<double >::lowest ();
483- }
484-
485- for (int corner = 0 ; corner < 8 ; corner++) {
486- Coord index;
487- index[0 ] = (corner & 1 ) ? size[0 ] - 1 : 0 ;
488- index[1 ] = (corner & 2 ) ? size[1 ] - 1 : 0 ;
489- index[2 ] = (corner & 4 ) ? size[2 ] - 1 : 0 ;
490-
491- Point3 point;
492- itk_image_->TransformIndexToPhysicalPoint (index, point);
493-
494- for (int dim = 0 ; dim < 3 ; dim++) {
495- minPt[dim] = std::min (minPt[dim], point[dim]);
496- maxPt[dim] = std::max (maxPt[dim], point[dim]);
497- }
498- }
499-
500- ImageType::DirectionType identityDir;
501- identityDir.SetIdentity ();
502-
503- ImageType::SizeType outputSize;
504- for (int i = 0 ; i < 3 ; i++) {
505- outputSize[i] = static_cast <size_t >(std::ceil ((maxPt[i] - minPt[i]) / spacing[i])) + 1 ;
506- }
507-
508- Vector3 spacingVec;
509- Dims dims;
510- for (int i = 0 ; i < 3 ; i++) {
511- spacingVec[i] = spacing[i];
512- dims[i] = outputSize[i];
513- }
514-
515- return resample (transform, minPt, dims, spacingVec, identityDir, interp);
516- }
517-
518436bool Image::compare (const Image& other, bool verifyall, double tolerance, double precision) const {
519437 if (tolerance > 1 || tolerance < 0 ) {
520438 throw std::invalid_argument (" tolerance value must be between 0 and 1 (inclusive)"
0 commit comments