View source for Coordinate System Conversion Between ITK and Slicer3

The coordinate systems used in ITK and Slicer3 are different. Internally in Slicer3 a transform is built that converts from an ijk pixel index into a RAS coordinate system. This is extracted from the image using the Direction Cosines from an ITK image with the GetDirection() method. A problem for the meshing programs (solid and surface based meshes) is that we need to extract the points in an RAS system coordinate system for the node locations. A reasonable simple way to handle this is to use the OrientedImage image type in ITK. However the Direction Cosines are for an LPS coordinate system and thus we need to convert the Direction Cosines to place place the image into a RAS coordinate system. Based on the code from the OrientImageFilter in ITK, a conversion of the Direction Cosines to RAS was designed that allows all of the coordinate transforms to be handled internally in the OrientedImage using the TransformIndexToPhysicalPoint() or   TransformContinuousIndexToPhysicalPoint() methods. 

 

  ImageType::DirectionType imageDir = image->GetDirection( );
  ImageType::PointType origin = image->GetOrigin( );
  
  int dominantAxisRL = itk::Function::Max3(imageDir[0][0],imageDir[1][0],imageDir[2][0]);
  int signRL = itk::Function::Sign(imageDir[dominantAxisRL][0]);
  int dominantAxisAP = itk::Function::Max3(imageDir[0][1],imageDir[1][1],imageDir[2][1]);
  int signAP = itk::Function::Sign(imageDir[dominantAxisAP][1]);
  int dominantAxisSI = itk::Function::Max3(imageDir[0][2],imageDir[1][2],imageDir[2][2]);
  int signSI = itk::Function::Sign(imageDir[dominantAxisSI][2]);
  
  ImageType::DirectionType DirectionToRAS;
  DirectionToRAS.SetIdentity( );
  
  if (signRL == 1)
    {
      DirectionToRAS[dominantAxisRL][dominantAxisRL] = -1.0; 
      origin[dominantAxisRL] *= -1.0;
    }
  if (signAP == -1)
    {
      DirectionToRAS[dominantAxisAP][dominantAxisAP] = -1.0; 
      origin[dominantAxisAP] *= -1.0;
    }
  if (signSI == 1)
    {
      DirectionToRAS[dominantAxisSI][dominantAxisSI] = -1.0; 
      origin[dominantAxisSI] *= -1.0;
    }
   
  imageDir *= DirectionToRAS;
 
  image->SetDirection( imageDir  );
  image->SetOrigin( origin );