convertToITK.cxx

//
// to ITK 3D
//
template <class T, unsigned int VImageDimension >
typename DeformationImageType::Pointer VectorImageUtils< T, VImageDimension >::convertToITK( const VectorImageType3D* im)
{

    /* DeformationImageType instead of ITKVectorImage<T,VImageDimension>::Type
     **
     ** We defined DeformationImageType like this in VectorImageUtils.h :
     ** typedef itk::Vector<T,VImageDimension> DeformationPixelType;
     ** typedef itk::Image<DeformationPixelType, VImageDimension> DeformationImageType;
     **
     */

  if ( VImageDimension != 3 )
    {
    throw std::runtime_error( "convertToITK3D only for 3D images." );
    }

  unsigned int szX = im->GetSizeX();
  unsigned int szY = im->GetSizeY();
  unsigned int szZ = im->GetSizeZ();
  unsigned int dim = im->GetDimension();

  // Initialize ITK image
  typename DeformationImageType::Pointer outImage;
  outImage = DeformationImageType::New();

  // Set up region
  typename DeformationImageType::IndexType start;
  start[0] = 0;
  start[1] = 0;
  start[2] = 0;
  //start[3] = 0;

  typename DeformationImageType::SizeType size;
  size[0] = szX;
  size[1] = szY;
  size[2] = szZ;
  //size[3] = dim;

  typename DeformationImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);

  // Set up the spacing
  typename DeformationImageType::SpacingType space;
  space[0] = im->GetSpacingX();
  space[1] = im->GetSpacingY();
  space[2] = im->GetSpacingZ();
  //space[3] = 1;
  outImage->SetSpacing(space);

  // Allocate region to image
  outImage->SetRegions(region);
  outImage->Allocate();

  // Copy in the data
  for (unsigned int z = 0; z < szZ; ++z)
    {
    for (unsigned int y = 0; y < szY; ++y)
      {
      for (unsigned int x = 0; x < szX; ++x)
        {
          typename DeformationImageType::IndexType px;
          DeformationPixelType pixel;
          px[0] = x;
          px[1] = y;
          px[2] = z;

        for (unsigned int d = 0; d < dim; ++d)
          {

            pixel[d] = im->GetValue(x,y,z,d);
            
          //typename DeformationImageType::Type::IndexType px;
          //px[0] = x;
          //px[1] = y;
          //px[2] = z;
          //px[3] = d;

          //outImage->SetPixel(px, im->GetValue(x,y,z,d));
          }
          
          outImage->SetPixel(px, pixel);
          
        }
      }
    }

  // Set origin and direction
  outImage->SetOrigin(im->GetOrigin());         //should not be ok (4D/3D)
  outImage->SetDirection(im->GetDirection());   //should not be ok (4D/3D)

  // return the result
  return outImage;

}