ENH: Add DWI scan padding

Now pads image on all sides

AutoWorkup/SEMTools/diffusion/gtract.py

@@ -714,12 +714,13 @@ class gtractResampleDWIInPlaceInputSpec(CommandLineInputSpec):
     inputVolume = File(desc="Required: input image is a 4D NRRD image.", exists=True, argstr="--inputVolume %s")
     inputTransform = File(desc="Required: transform file derived from rigid registration of b0 image to reference structural image.", exists=True, argstr="--inputTransform %s")
     debugLevel = traits.Int(desc="Display debug messages, and produce debug intermediate results.  0=OFF, 1=Minimal, 10=Maximum debugging.", argstr="--debugLevel %d")
-    outputVolume = traits.Either(traits.Bool, File(), hash_files=False, desc="Required: output image (NRRD file) that has been transformed into the space of the structural image.", argstr="--outputVolume %s")
+    imageOutputSize = InputMultiPath(traits.Int, desc="The voxel lattice for the output image, padding is added if necessary. NOTE: if 0,0,0, then the inputVolume size is used.", sep=",", argstr="--imageOutputSize %s")
+    outputVolume = traits.Either(traits.Bool, File(), hash_files=False, desc="Required: output image (NRRD file) that has been rigidly transformed into the space of the structural image and padded if image padding was changed from 0,0,0 default.", argstr="--outputVolume %s")
     numberOfThreads = traits.Int(desc="Explicitly specify the maximum number of threads to use.", argstr="--numberOfThreads %d")
 
 
 class gtractResampleDWIInPlaceOutputSpec(TraitedSpec):
-    outputVolume = File(desc="Required: output image (NRRD file) that has been transformed into the space of the structural image.", exists=True)
+    outputVolume = File(desc="Required: output image (NRRD file) that has been rigidly transformed into the space of the structural image and padded if image padding was changed from 0,0,0 default.", exists=True)
 
 
 class gtractResampleDWIInPlace(SEMLikeCommandLine):
@@ -735,7 +736,7 @@ class gtractResampleDWIInPlace(SEMLikeCommandLine):
 
 license: http://mri.radiology.uiowa.edu/copyright/GTRACT-Copyright.txt
 
-contributor: This tool was developed by Vincent Magnotta and Greg Harris.
+contributor: This tool was developed by Vincent Magnotta, Greg Harris, Hans Johnson, and Joy Matsui.
 
 acknowledgements: Funding for this version of the GTRACT program was provided by NIH/NINDS R01NS050568-01A2S1
 

AutoWorkup/SEMTools/registration/brainsfit.py

@@ -23,6 +23,7 @@ class BRAINSFitInputSpec(CommandLineInputSpec):
     useComposite = traits.Bool(desc="Perform a Composite registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useComposite ")
     numberOfSamples = traits.Int(desc="The number of voxels sampled for mutual information computation.  Increase this for a slower, more careful fit.  You can also limit the sampling focus with ROI masks and ROIAUTO mask generation.", argstr="--numberOfSamples %d")
     splineGridSize = InputMultiPath(traits.Int, desc="The number of subdivisions of the BSpline Grid to be centered on the image space (in voxels).  Each dimension must have at least 3 subdivisions for the BSpline to be correctly computed. ", sep=",", argstr="--splineGridSize %s")
+    useROIBSpline = traits.Bool(desc="When set, will use the bounding box of the input ROIs to define spline grid support region. Otherwise will use the whole image.", argstr="--useROIBSpline ")
     numberOfIterations = InputMultiPath(traits.Int, desc="The maximum number of iterations to try before failing to converge.  Use an explicit limit like 500 or 1000 to manage risk of divergence", sep=",", argstr="--numberOfIterations %s")
     maskProcessingMode = traits.Enum("NOMASK", "ROIAUTO", "ROI", desc="What mode to use for using the masks.  If ROIAUTO is choosen, then the mask is implicitly defined using a otsu forground and hole filling algorithm. The Region Of Interest mode (choose ROI) uses the masks to define what parts of the image should be used for computing the transform.", argstr="--maskProcessingMode %s")
     fixedBinaryVolume = File(desc="Fixed Image binary mask volume, ONLY FOR MANUAL ROI mode. (Values are zero and non-zero.)", exists=True, argstr="--fixedBinaryVolume %s")

AutoWorkup/SEMTools/registration/brainsfitez.py

@@ -19,9 +19,11 @@ class BRAINSFitEZInputSpec(CommandLineInputSpec):
     useScaleVersor3D = traits.Bool(desc="Perform a ScaleVersor3D registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useScaleVersor3D ")
     useScaleSkewVersor3D = traits.Bool(desc="Perform a ScaleSkewVersor3D registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useScaleSkewVersor3D ")
     useBSpline = traits.Bool(desc="Perform a BSpline registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useBSpline ")
+    useSyN = traits.Bool(desc="Perform a SyN registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useSyN ")
     useComposite = traits.Bool(desc="Perform a Composite registration as part of the sequential registration steps.  This family of options superceeds the use of transformType if any of them are set.", argstr="--useComposite ")
     numberOfSamples = traits.Int(desc="The number of voxels sampled for mutual information computation.  Increase this for a slower, more careful fit.  You can also limit the sampling focus with ROI masks and ROIAUTO mask generation.", argstr="--numberOfSamples %d")
     splineGridSize = InputMultiPath(traits.Int, desc="The number of subdivisions of the BSpline Grid to be centered on the image space (in voxels).  Each dimension must have at least 3 subdivisions for the BSpline to be correctly computed. ", sep=",", argstr="--splineGridSize %s")
+    useROIBSpline = traits.Bool(desc="When set, will use the bounding box of the input ROIs to define spline grid support region. Otherwise will use the whole image.", argstr="--useROIBSpline ")
     numberOfIterations = InputMultiPath(traits.Int, desc="The maximum number of iterations to try before failing to converge.  Use an explicit limit like 500 or 1000 to manage risk of divergence", sep=",", argstr="--numberOfIterations %s")
     maskProcessingMode = traits.Enum("NOMASK", "ROIAUTO", "ROI", desc="What mode to use for using the masks.  If ROIAUTO is choosen, then the mask is implicitly defined using a otsu forground and hole filling algorithm. The Region Of Interest mode (choose ROI) uses the masks to define what parts of the image should be used for computing the transform.", argstr="--maskProcessingMode %s")
     fixedBinaryVolume = File(desc="Fixed Image binary mask volume, ONLY FOR MANUAL ROI mode. (Values are zero and non-zero.)", exists=True, argstr="--fixedBinaryVolume %s")

GTRACT/Cmdline/gtractResampleDWIInPlace.cxx

@@ -96,6 +96,8 @@ int main(int argc, char *argv[])
     std::cout << "Input Transform: " <<  inputTransform << std::endl;
     std::cout << "Output Image: " <<  outputVolume << std::endl;
     std::cout << "Debug Level: " << debugLevel << std::endl;
+    std::cout << "Image Output Size: "
+      << imageOutputSize[0] << "," << imageOutputSize[1] << "," << imageOutputSize[2] << std::endl;
     std::cout << "=====================================================" << std::endl;
     }
 
@@ -120,7 +122,7 @@ int main(int argc, char *argv[])
   typedef signed short                        PixelType;
   typedef itk::VectorImage<PixelType, 3>      NrrdImageType;
   typedef itk::VersorRigid3DTransform<double> RigidTransformType;
-
+  typedef itk::Image<PixelType, 3>            InputIndexImageType;
   typedef itk::ImageFileReader<NrrdImageType,
                                itk::DefaultConvertPixelTraits<PixelType> > FileReaderType;
   FileReaderType::Pointer imageReader = FileReaderType::New();
@@ -162,8 +164,9 @@ int main(int argc, char *argv[])
 
   // Get measurement frame and its inverse from DWI scan
   std::vector< std::vector< double > > msrFrame;
-  itk::ExposeMetaData< std::vector< std::vector< double > > >( inputMetaDataDictionary, "NRRD_measurement frame", msrFrame );
-
+  itk::ExposeMetaData< std::vector< std::vector< double > > >( inputMetaDataDictionary,
+    "NRRD_measurement frame", msrFrame );
+
 
   vnl_matrix_fixed< double, 3, 3 > DWIMeasurementFrame;
   for( unsigned int i = 0; i < 3; i++ )
@@ -173,16 +176,15 @@ int main(int argc, char *argv[])
       DWIMeasurementFrame[ i ][ j ] = msrFrame[ i ][ j ];
     }
   }
-  std::cout << "DWI measurement frame (DWIMeasurementFrame): " << DWIMeasurementFrame << std::endl;
   vnl_matrix_fixed< double, 3, 3 > DWIInverseMeasurementFrame = vnl_inverse( DWIMeasurementFrame );
-  std::cout << "DWI inverse measurement frame (DWIInverseMeasurementFrame): " << DWIInverseMeasurementFrame << std::endl; 
-
-
+
+
   // Resample DWI in place
   resampleImage = SetVectorImageRigidTransformInPlace<NrrdImageType>(rigidTransform.GetPointer(), resampleImage);
 
-
   std::cout << "Rigid transform matrix: " << rigidTransform->GetMatrix().GetVnlMatrix() << std::endl;
+
+
   // Rotate gradient vectors by rigid transform and inverse measurement frame
   for( unsigned int i = 0; i < resampleImage->GetVectorLength(); i++ )
     {
@@ -198,11 +200,8 @@ int main(int argc, char *argv[])
       NrrdValue.c_str(), " %lf %lf %lf", &curGradientDirection[0], &curGradientDirection[1], &curGradientDirection[2]);
 
     // Rotate the diffusion gradient with rigid transform and inverse measurement frame
-    std::cout << "Current gradient direction (before rotation): " << curGradientDirection << std::endl;
     RigidTransformType::Pointer inverseRigidTransform = RigidTransformType::New();
-    const bool invertWorked = rigidTransform->GetInverse( inverseRigidTransform );
     curGradientDirection = inverseRigidTransform->GetMatrix().GetVnlMatrix() * DWIInverseMeasurementFrame * curGradientDirection;
-    std::cout << "Current gradient direction (after rotation): " << curGradientDirection << std::endl;
 
     // Updated the Image MetaData Dictionary with Updated Gradient Information
     //sprintf(tmpStr, " %18.15lf %18.15lf %18.15lf", curGradientDirection[0], curGradientDirection[1],
@@ -220,7 +219,6 @@ int main(int argc, char *argv[])
   // Set DWI measurement frame to identity by multiplying by its inverse
   // Update Image MetaData Dictionary with new measurement frame
   vnl_matrix_fixed< double, 3, 3 > newMeasurementFrame = DWIInverseMeasurementFrame * DWIMeasurementFrame;
-  std::cout << "New measurement frame (newMeasurementFrame): " << newMeasurementFrame << std::endl;
   std::vector< std::vector< double > > newMf(3);
   for( unsigned int i = 0; i < 3; i++ )
   {
@@ -230,15 +228,119 @@ int main(int argc, char *argv[])
       newMf[ i ][ j ] = newMeasurementFrame[ i ][ j ];
     }
   }
-  itk::EncapsulateMetaData< std::vector< std::vector< double > > >( resampleImage->GetMetaDataDictionary(), "NRRD_measurement frame", newMf ); 
+  itk::EncapsulateMetaData< std::vector< std::vector< double > > >( resampleImage->GetMetaDataDictionary(), "NRRD_measurement frame", newMf );
+
+
+  // Pad image
+  const NrrdImageType::RegionType inputRegion = resampleImage->GetLargestPossibleRegion();
+  const NrrdImageType::SizeType inputSize = inputRegion.GetSize();
+  const NrrdImageType::PointType inputOrigin = resampleImage->GetOrigin();
+  const NrrdImageType::DirectionType inputDirection = resampleImage->GetDirection();
+  const NrrdImageType::SpacingType inputSpacing = resampleImage->GetSpacing();
+
+  NrrdImageType::SizeType newSize;
+  std::vector<size_t> imagePadding(6,0);
+  for( int qq =0; qq < 3 ; ++qq )
+    {
+    if ( ( imageOutputSize[qq] > 0 ) && ( static_cast<itk::SizeValueType>( imageOutputSize[qq] )  > inputSize[qq] ) )
+      {
+      const size_t sizeDiff = imageOutputSize[qq] - inputSize[qq];
+      const size_t halfPadding = sizeDiff/2;
+      const size_t isOddDiff = sizeDiff%2;
+      imagePadding[qq] =  halfPadding  + isOddDiff;
+      imagePadding[qq+3] = halfPadding;
+      }
+     newSize[qq] = inputSize[qq] + imagePadding[qq] + imagePadding[qq+3];
+    }
+
+  vnl_matrix_fixed< double, 3, 3 > inputDirectionMatrix;
+  vnl_matrix_fixed< double, 3, 3 > inputSpacingMatrix;
+
+  for( unsigned int i = 0; i < 3; i++ )
+  {
+    for ( unsigned int j = 0; j < 3; j++ )
+    {
+      inputDirectionMatrix[ i ][ j ] = inputDirection[ i ][ j ];
+      if( i == j )
+      {
+        inputSpacingMatrix[ i ][ j ] = inputSpacing[ i ];
+      }
+      else
+      {
+        inputSpacingMatrix[ i ][ j ] = 0.0;
+      }
+    }
+  }
+
+  vnl_matrix_fixed< double, 3, 3 > spaceDirections = inputDirectionMatrix * inputSpacingMatrix;
+  vnl_matrix_fixed< double, 3, 4 > newMatrix;
+  for( unsigned int i = 0; i < 3; i++ )
+  {
+    for( unsigned int j = 0; j < 4; j++ )
+    {
+      if( j == 3 )
+      {
+        newMatrix[ i ][ j ] = inputOrigin[ i ];
+      }
+      else
+      {
+        newMatrix[ i ][ j ] = spaceDirections[ i ][ j ];
+      }
+    }
+  }
+
+  vnl_matrix_fixed< double, 4, 1 > voxelShift;
+  voxelShift[ 0 ][ 0 ] = -1.0 * ( double )( imagePadding[ 0 ] );
+  voxelShift[ 1 ][ 0 ] = -1.0 * ( double )( imagePadding[ 1 ] );
+  voxelShift[ 2 ][ 0 ] = -1.0 * ( double )( imagePadding[ 2 ] );
+  voxelShift[ 3 ][ 0 ] = 1.0;
+
+  vnl_matrix_fixed< double, 3, 1 > newOriginMatrix = newMatrix * voxelShift;
+
+  NrrdImageType::PointType newOrigin;
+  for( unsigned int i = 0; i < 3; i++ )
+  {
+    newOrigin[ i ] = newOriginMatrix[ i ][ 0 ];
+  }
+
+  std::cout << "Input DWI Image Origin: ( " << inputOrigin[ 0 ] << ", " << inputOrigin[ 1 ] << ", " << inputOrigin[ 2 ] << " )" << std::endl;
+  std::cout << "Input DWI Image Size: " << inputSize[ 0 ] << " " << inputSize[ 1 ] << " " << inputSize[ 2 ] << std::endl;
+  std::cout << " " << std::endl;
+  std::cout << "Output DWI Image Origin: ( " << newOrigin[ 0 ] << ", " << newOrigin[ 1 ] << ", " << newOrigin[ 2 ] << " )" << std::endl;
+  std::cout << "Output DWI Image Size: " << newSize[ 0 ] << " " << newSize[ 1 ] << " " << newSize[ 2 ] << std::endl;
+
+
+  NrrdImageType::Pointer paddedImage = NrrdImageType::New();
+  paddedImage->CopyInformation(resampleImage);
+  paddedImage->SetVectorLength( resampleImage->GetVectorLength() );
+  paddedImage->SetMetaDataDictionary( resampleImage->GetMetaDataDictionary() );
+  paddedImage->SetRegions( newSize );
+  paddedImage->SetOrigin( newOrigin );
+  paddedImage->Allocate();
+
+  typedef itk::ImageRegionIterator< NrrdImageType > IteratorType;
+  IteratorType InIt( resampleImage, resampleImage->GetRequestedRegion() );
+
+  for( InIt.GoToBegin(); !InIt.IsAtEnd(); ++InIt )
+  {
+    NrrdImageType::IndexType InIndex = InIt.GetIndex();
+    NrrdImageType::IndexType OutIndex;
+    OutIndex[ 0 ] = InIndex[ 0 ] + imagePadding[ 0 ];
+    OutIndex[ 1 ] = InIndex[ 1 ] + imagePadding[ 1 ];
+    OutIndex[ 2 ] = InIndex[ 2 ] + imagePadding[ 2 ];
+
+    NrrdImageType::PixelType InImagePixel = resampleImage->GetPixel( InIndex );
+    paddedImage->SetPixel( OutIndex, InImagePixel );
+    paddedImage->Update();
+  }
+
 
-
   // Write out resampled in place DWI
   typedef itk::ImageFileWriter<NrrdImageType> WriterType;
   WriterType::Pointer nrrdWriter = WriterType::New();
   nrrdWriter->UseCompressionOn();
   nrrdWriter->UseInputMetaDataDictionaryOn();
-  nrrdWriter->SetInput( resampleImage );
+  nrrdWriter->SetInput( paddedImage );
   nrrdWriter->SetFileName( outputVolume );
   try
     {

GTRACT/Cmdline/gtractResampleDWIInPlace.xml

@@ -8,7 +8,7 @@
   <version>4.0.0</version>
   <documentation-url>http://wiki.slicer.org/slicerWiki/index.php/Modules:GTRACT</documentation-url>
   <license>http://mri.radiology.uiowa.edu/copyright/GTRACT-Copyright.txt</license>
-  <contributor>This tool was developed by Vincent Magnotta and Greg Harris.</contributor>
+  <contributor>This tool was developed by Vincent Magnotta, Greg Harris, Hans Johnson, and Joy Matsui.</contributor>
 
 
   <parameters>
@@ -41,14 +41,27 @@
 
   </parameters>
 
+  <parameters advanced="true">
+    <label>Main Parameters</label>
+
+    <integer-vector>
+      <name>imageOutputSize</name>
+      <longflag>imageOutputSize</longflag>
+      <label>Image Output Padding</label>
+      <description>The voxel lattice for the output image, padding is added if necessary. NOTE: if 0,0,0, then the inputVolume size is used.</description>
+      <default>0,0,0</default>
+    </integer-vector>
+
+  </parameters>
+
   <parameters>
     <label>Output Files</label>
     <description></description>
 
     <image type="diffusion-weighted">
       <name>outputVolume</name>
       <longflag>outputVolume</longflag>
-      <description>Required: output image (NRRD file) that has been transformed into the space of the structural image.</description>
+      <description>Required: output image (NRRD file) that has been rigidly transformed into the space of the structural image and padded if image padding was changed from 0,0,0 default.</description>
       <label>Output File</label>
       <channel>output</channel>
     </image>